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ABSTRACT 


An  experimental  computer  program  is  described,  which  calculates 
Earth  Satellite  ephemerides,  corrects  orbit  elements  and  evaluates  the 
effects  of  various  terms  of  the  bulge  perturbation  theory.  Perturbations 
by  solar  radiation  pressure  and  atmospheric  drag  are  also  represented. 

The  differential  correction  employs  a  weighted  least-squares  reduction. 
Formulation,  flow  charts,  input  formats  and  sample  cases  are  given. 
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SECTION  1 

INTRODUCTION 


The  Experimental  General  Perturbations  Differential  Correction 
Program  calculates  the  ephemeris  of  an  Earth  satellite  by  a  General 
Perturbations  technique  and  improves  the  orbital  elements  by  a  differential 
correction  process  using  satellite  observations.  The  ephemeris  calculation 
includes  the  analytical  expressions  of  the  perturbations  caused  by  the 
asphericity  of  the  Earth  and  of  the  effects  due  to  direct  solar  radiation 
pressure  on  a  close-Earth  satellite.  Frictional  effects  due  to  the  Earth's 
atmosphere  are  determined  empirically  from  observational  data. 

The  program  is  equipped  to  compare  the  complete  first-order 
asphericity  theory  with  simplified  theories,  in  which  selected  terms  are 
omitted.  The  comparison  is  made  in  the  magnitude  and  in  the  radial,  trans¬ 
verse  and  orthogonal  components  of  the  displacement  from  the  position 
obtained  with  the  complete  theory. 

The  effects  of  direct  solar  radiation  on  the  orbit  of  an  Earth 
satellite  are  introduced  through  the  perturbations  in  the  orbital  parameters 


and 


e 

cos 

7T 

e 

sin 

7 T 

L,  mean  longitude 

sin  i  sin 
sin  i  cos 


-1- 


doT^J^tO^SoTTl/ia^iy-, 

AERONUTRONIC  DIVISION 


where 


and 


7r  =  +  0)  (see  also  Figure  3) 

ft  is  the  longitude  of  the  ascending  node 
co  is  the  argument  of  perigee 
e  is  the  eccentricity  of  orbit 

i  is  the  inclination  of  the  orbit— plane  to  the  equator-plane. 


This  set  of  elements  defines  nearly -circular  orbits  without  the  difficulties 
caused  by  the  singularities  of  the  classical  elements  at  zero  eccentricity. 

Only  the  long-period  changes  of  the  parameters  are  retained  from  the  develop¬ 
ment.*  The  effects  on  these  long-period  terms  of  the  satellite  being  eclipsed 
by  the  Earth  in  part  of  its  orbit  are  also  included. 

The  perturbation  theory  for  the  asphericity  of  the  Earth  contains 

first-order  short-period  terms  and  second-order  secular  and  long-period  terms. 

The  second  zonal  harmonic  coefficient,  J^,  of  the  geopotential  function  is 

considered  to  be  of  the  first  order.  The  third  and  fourth  zonal  harmonic 

o 

coefficients,  and  J^,  respectively,  as  well  as  J  are  of  the  second  order. 

For  a  more  detailed  analysis,  reference  is  rade  to  Aeronutronic  Report  S-981. 

In  addition,  the  works  of  Brouwer**  and  Kozai***  are  useful  for  comparison  of 
the  results .  The  second-order  secular  and  long-period  terms  in  the  mean 
anomaly,  M,  were  adopted  from  Brouwer's  analysis. 

The  expressions  presented  in  the  formulation  of  these  bulge  perturba¬ 
tions  are  free  of  low-eccentricity  singularities  due  to  the  choice  of  the 
parameters  e  sin  CO,  e  cos  co  and  the  argument  of  latitude,  u,  instead  of  the 
corresponding  classical  elements.  The  short-period  terms  in  the  radial 
distance,  r,  and  the  argument  of  latitude,  u,  are  used  rather  than  the  related 
expressions  for  the  semi -axis  major,  a,  the  eccentricity,  e,  the  argument  of 
perigee,  co  ,  and  the  mean  anomaly,  M,  which  contain  ^  .  Long -period  terms  are 
calculated  for  the  elements  e  cos  CO,  e  sin  CO,  and  L  since  the  classical  parameters 
also  have  terms  containing  J^/e  as  a  coefficient. 


*  Koskela,  Paul,  1961,  "Orbital  Effects  of  Solar  Radiation  Pressure 
on  an  Earth  Satellite,"  Aeronutronic  Publication  No.  U-1357. 

**  Brouwer,  Dirk,  1959,  Astronomical  Journal,  64.  378. 

***  Kozai,  Y,  1959,  loc .  cit,  367. 
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Several  of  the  bulge  perturbation  terms  contain  the  factor 


which  approaches  %  as  a  limit  at  zero  eccentricity.  It  is  approximated  by  a 
rapidly  converging  series  expansion  when  e  is  small,  as  are  two  similar  co¬ 
efficients  . 

2 

The  expression  (4-5  sin  i)  appears  in  the  denominator  of  many  of  the 
long-period  terms  shown  in  the  formulation.  This  quantity  produces  a  singu¬ 
larity  at  the  critical  inclination  angle  of  63 .4349 .. .degrees .  In  this 
program  the  terms  containing  (4-5  sin^i)  are  automatically  set  equal  to  zero 
when  the  inclination  approaches  the  critical  value.* 

The  perturbation  theory  is  based  on  the  Gaussian  form  of  the  expressions 
for  the  variations  of  the  parameters;  that  is,  three  mutually  perpendicular 
components  of  the  perturbing  acceleration  are  used.  The  expressions  are  de¬ 
veloped  in  terms  of  the  true  anomaly  and  thus  have  no  remainder  terms  in  powers 
of  the  eccentricity. 

Section  2  of  this  report  lists  the  equations  necessary  for  ephemeris 
calculation.  The  equations  are  presented  in  a  form  which  allows  clear  identi¬ 
fication  of  each  bulge  term  for  the  experimentation.  The  use  of  low-eccentricity 
orbital  parameters,  the  inclusion  of  the  effects  of  the  Earth's  bulge  on  the 
satellite's  velocity,  and  the  separation  of  the  terms  by  size  with  regard  to 
the  order  of  the  harmonic  coefficients  and  the  power  of  the  eccentricity  account 
for  the  large  number  of  terms . 


Section  3  describes  the  computer  program  which  has  resulted  from 
this  project.  The  ephemeris  calculation  consists  of  determining  the  pertur¬ 
bative  changes  in  the  selected  orbital  parameters  or  coordinates  and  calculating 
the  geocentric  position  and  velocity  vectors  at  desired  times.  The  differential 
correction  process  obtains  the  correction  to  any  or  all  of  the  six  orbital 

elements;  U  ,  the  mean  argument  of  latitude,  a  „  =  e  cos  CD  ,  a  „  =  e  sin  ffl  , 
o’  °  ’  xN  yN 

n  ,  the  orbital  mean  motion,  i  and  Q  ;  and  the  correction  to  the  drag  parameters 

c"  and  d,  if  desired.  The  program  has  the  capacity  to  weight  the  observations 

according  to  the  given  value  of  the  standard  deviation  of  each  observation,  be 

it  range,  azimuth,  elevation  angle,  right  ascension,  declination  or  range-rate. 

It  may  also  consider  the  weights  assigned  to  an  observing  station  for  the 

quality  of  its  observational  data. 


*  Aeronutronic  Report  U-1672,  "Numerical  Values  of  General  Perturbations  Terms 
in  Orbital  Parameters,"  by  J .  L.  Arsenault  and  its  appendix  explain  the  justi¬ 
fication  for  this.  Basically,  as  the  rate  of  change  in  the  argument  of  perigee 
approaches  zero,  these  particular  long-period  terms  become  of  very  long  period 
and  act  like  secular  terms.  The  mean  value  of  the  terms  can  be  successfully 
combined  with  the  mean  value  of  the  relevant  parameter. 
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SECTION  2 

FORMULATION 


Because  the  treatment  of  perturbations  due  to  direct  solar  radiation 
in  this  program  is  somewhat  novel,  the  formulation  is  detailed  here.  The 
formulation  of  all  the  ephemeris  calculation,  including  the  other  perturbations 
and  of  the  differential  correction  are  then  given. 


2 . 1  DIRECT  SOIAR  RADIATION  PERTURBATIONS 

The  magnitude  of  the  force  acting  on  a  satellite  due  to  direct  solar 
radiation  is 

F  =  I  F  |  =  -y  v  P  A 

where 

A  is  the  effective  cross-sectional  area  of  the  satellite 
to  radiation  pressure 

P  is  the  solar  radiation  pressure  in  the  vicinity  of  the 
Earth,  assumed  constant 

7  is  a  factor  depending  on  the  reflecting  characteristics 
of  the  satellite's  surface 

V  is  an  eclipse  factor  to  be  considered  in  detail  later 
(Sec t ion  2.1.3) 

It  is  customary  to  assume  that  the  direction  and  magnitude  of  this  force  is 
constant  along  the  entire  orbit,  except  in  the  Earth's  shadow,  where  it 
does  not  exist. 
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2.1.1  Sun's  True  Longitude 


The  Sun's  true  longitude,  at  time  t  is  computed  from 

l  «  (L  )  +  n  (t-t  )  +  2  e  sin  M  +  y  e?  sin  2  M  +  . . . (in  radians) 
e  •  o  •  o  a  o4o  o 


where 


(I/)  is  the  Sun's  mean  longitude  at  some  epoch,  tQ 

nQ  =  0?9856/day,  the  Sun's  mean  daily  motion 

e^  *  0.016726,  the  eccentricity  of  the  Earth's  orbit 

M  =  (L)  +  n  (t-t)  -  K 

o  o  o  o  o  o 

n  =  282°.  253 

o 


2.1.2  Perturbations  in  the  Elements  Due  to  Direct  Solar  Radiation 

Compute  the  perturbations  in  the  elements  a,  (e  cos  tv ) , 

(e  sin  tt) ,  W  ,  W  ,  and  L  due  to  direct  solar  radiation.  The  required 
expressions  are  derived  in  Aeronutronic  Report  U-1357.  However,  in  the 

d  fl  d  CD 

results  which  follow,  -rr~  and  — — -  have  been  replaced  by  their  variations 

dE  aE 

with  respect  to  the  modified  time  variable, T ,  by  means  of 


=  -s.  (I 

dE  n 


^  oi 
dE  n 

/  3 /2 

The  mean  motion,  n,  has  been  replaced  by  n  =  kg  *» /a  radians/min. 

The  a^  coefficients  are  listed  in  Table  I.  The  semi -major  axis 


a,  is  in  units  of  Earth  radii;  pt  =  1  radlus) —  .  t^e  satellite's 

(k  min) 
e 
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mass,m,  is  In  grams;  the  magnitude  of  the  force  due  to  direct  solar 
radiation  is  F  =  V'yp  A;  the  eclipse  factor  o  and  the  reflectivity  y 
are  dimensionless;  the  satellite's  effective  cross-sectional  area.  A, 

2  -8  n  radii 

is  in  cm  ;  and  the  solar  radiation  pressure  is  T  *  4.592  x  10  — — r—  _ 

e  cm  (kg  min) 

(a)  In  the  semi-major  axis  there  are  only  short-period  perturbations 
due  to  direct  solar  radiation.  Therefore,  let 

A  a  *  0 


Define  an  unperturbed  semi-axis  major: 

3  J0a  2  , - 

a  »  a  [1+  - =r”  1-e 

uo  o  1  „  2  o 


a  =«a  [1+  - 

UP  °  2  p 

* 


*/  1-e  2  (1  -  3/2  sin  2i  )] 

o  o 


(b)  Calculate  A  (e  cos  7r) : 

A  (e  cos  it)  =  -  T?  f  &  _  (1-e2)  J  4  sin  i  sin  e 

16  UP  °  sTTTm  [ 

[  -  a  cos  (ft  +  i  )  +  a  cos  (ft  -  l  )] 


-  2  (1  -  cos  iQ)  (1  ■  cose)  cos  (2  ft+  /^) 


+  a14  (1  +  cos  e)  cos  (2ft-  i#)  ] 


-4  a..  (1  +  cos  i  )  (  cos  e  )  cos  / 


sin  i  F  ( 

_ o  o  / 

+  cos  iQ)  •/'jjm  ( 


2  cos  i  sin  € 
o 


r  -  a  cos  (ft+  /  ) 
1  11  o 


+  a12  cos  (ft  -  i0)  +  a16  cos  (ft+  2cu  +  /@) 


-  a17  cos  (ft+  2o)-  /e)  ] 
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+  sin  iQ  (1  +  cos  e  )  [  -  a14  cos  (2  ft  -  £  0) 

-  alg  cos  (2<jd  +  £J  +  a2Q  cos  (2  ft  +  2  0)  -  i  ] 

+  sin  i  (1  -  cos  €  )  [  -  a  cos  (2  ft  +  i  ) 

o  i  J  • 

-  a19  cos  (2  cd  -  i@)  +  «21  cos  (2  ft  +  2cjd  +  IJ  ] 


+  2  a,„  sin  i  cos  e  cos  £ 
15  o  a 


(c)  Calculate  A  (e  sin  7r  ): 


I 


A  (e  sinTT)  -  -  ^  ^d-e^ 


I 


4  sin  i  sin  £ 
o 


[  -a11  sin  (ft  +  /#)  +  a12  sin  (ft  -  i#)  ] 

-  2  (1  -  cos  i  )  [a  (1  -  cose)  sin  (2ft+i  ) 

O  IJ  • 


+  a14  (1  +  cos  e)  sin  (2ft  ~i#)] 


4  a^  (1  +  cos  i^)  sin  £ 


sT 


•) 


3  au„  eo  8in  1 


3  up  u  o  F«  r_ 

_  f  .  -  - <  2  cos  i 

*  1  -  e  (1  +  cos  i  )  * 


sin  £ 


[-  au  sin  (ft+  £J 


+  a12  sin  (  ft-  i@)  +a16  sin  ( ft+  2cd+ 


-  a17  sin  (ft+  2cu  -  £  )  ] 


+  sin  i  (1  +  cos  e)[  -  a.,  sin  (2ft-  £  ) 
o  1  14  » 


-  a1Q  sin  (2cd  +  i#)  +  a2Q  sin  (2 ft+  2CD-  £  J)  ] 
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+  sin  iQ  (I  -  cos  e)  [  -  a^  sin  (2(1+  i#) 

-  a19  sin  (20J -  tj  +  a21  sin  (20+200  +  /#)  ] 

+  2  a,,  sin  i  sin  i  l 
15  o  •  / 

(d)  Calculate  AW  : 

x 

^a  e  F  f 

A  WX  “  16  /7?"7  \  2  C0S  io  8in  6(1  +  cos  V 

o 

[  -  ax  sin  (fl+6D  +  £0)  +  *1  sin  (ft+o>-  /0)  ] 

+  2  cos  i  sin  e  (1  -  cos  i  )  [a_  sin  ((1-CD+  l  ) 
o  o  3  e 

-  a4  sin  (fl  -  CD-  ] 

+  2  a  sin  i  (cos  i  +  cos  e)  sin  (a>  + 1  ) 

5  o  o  e 

+  2  a,  sin  i  (cos  i  -  cos  e)  sin  (oj  ~  i  ) 
o  o  o  o 

+  a7  sin  iQ  (1  -  cos  iQ)  (1  +  cos  e)  sin  (2  ()- a)'  i  #) 

+  a  sin  i  (1  -  cos  i  )  (1  -  cos  e)  sin  (2  ft-o>+  /  ) 
o  o  o  o 

-  a  sin  i  (1  +  cos  i  )  (1  +  cose)  sin  (2fl+C0-  l  ) 

7  0  0  o 

-a1Q  sin  iQ  (1  +  cos  iQ)  (1  -  cos  e)  sin  (2fl +co+  j 

(e)  Calculate  AW  : 

y 

^  </au  eQ  F  J 

AWy  ■  r-T2  \  2  C0S  i0  Sin  e(l  +  COS  iQ)  [  •1  COS  (fl+CD+/0) 

o  ^ 

-  a2  cos  (£l  +  0>-  /#)  ] 
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(f)  Calculate  AL: 


3  >7 - 

-  —  a..  e 


+  2  coa  i  sin  e(l  -  cos  i  )  [-a  cob  (ft -£0  +1  ) 

O  O  J  • 

+  a^  cos  (ft-u>-  i#)] 

-  2  a„  sin  i  (1  +cos  i  cos  e)  cos  (o>+i  ) 

5  o  o  • 

-  2  a,  sin  i  (1  -  cos  i  cos  e)  cos  (o>-  /  ) 

boo  • 

-  a,  sin  i  (1  -  cos  i  )  (1  +  cose)  cos  (2ft-a>-  i  ) 

-  aa  sin  i  (1  -  cos  i  )  (1  -  cos  e)  cos  (2fl -co+  l  ) 

8  o  o  9 

+  a9  sin  iQ  (1  +  cos  iQ)  (1  +  cos  €)  cos  (2ft+0>-  IJ 

-> 

+  a,rt  sin  i  (1  +  cos  i  )  (1  -  cos  e)  cos  (2ft+0)  +  /  )  ? 
Woo  9  j 


o 


u  eQ  1 - J- ■  M  —  <  a  (1-tcos  iQ)  (1-cos  e)  sino(  ft-KO+/#) 

P  l  2(l+%/l-e^)jV7i’in 


+  a„  (1  +  cos  i  )  (1+  cose)  sin  (ft  +  co-  i  ) 

4  o  " 

+  a3  (1  -  cos  iQ)  (1  -  co3  e)  sin  (ft-co  +  i#) 

+  a.  (1  -  cos  i  )  (1+  cose)  sin  (ft-o>-  i  ) 

4  0  • 

1 

+  2  sin  iQ  sine  [  -  a^  sin  <  0)+  )  +  «6  sin  (  co-i^)  ]  j- 

_  V”u  eQ  sin  i  F  f 

+  r  - 2 -  —  <  sin  i  (1  -  cose)  [-a  sin  (ft  +  CD+Z  ) 

*'1l-e02  (1  +  cos  i  )  •/>«[  ° 

+  a^  sin  (ft  -<u  +  /.)] 

+  sin  i  (1  +  cos  e  )  [a,  sin  (ft-CD-  /  )  *  a„  sin  (ft+co-/  )] 

+  2  cos  i  sine  [  "  aj  sin  /0)  +  a$  ain  (  J  ^ 
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TABLE  I 

The  Coefficients  a 


-1 


vl 


radians/k  min 
e 


a,  -(rf+co'+n) 
1  0 

a2  *  (  ft'  +  0DX  -  n#) 

a  =  (  ftv  -  cox  +  n  ) 

J  © 

a,  -  (  ft'-  0)'-  n  ) 

4  © 

a  -  (  0)'  +  n  )  1 

5  © 

a.  *  (  0)x  -  n  )  * 

6  o 

a,  -  (2f£-0)x-  n  )’ 

aa  ■  (2^-0J+  n  ) 

O  0 

an  -  (2(t+CD-  n  )” 

a1Q  -  (2d+0$  +  n#) 

all  "  («N;fn0)'1 


al3 

al4 

a15 

a16 

a17 

al8 

a19 

a20 


\  _i 

(  n  -  n  ) 

0 

(2it+  no)_1 

(2d-  n  )’1 
0 

(tf+  2aJ  +  n  )-1 
0 

(d+  2a>-  n  )_1 

0 

(20)+  n^)"1 

(2  CD  -  no)_1 

(2  d+  2  oi  -  n  ) 
o 

(2d+  2CU+  n  )' 
o 


Note:  n 

e 


1.6064  x  10 


radians 

k  ^  min 
e 
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2.1.3  Eclipse  Factor 

When  a  satellite  Is  continuously  exposed  to  sunlight,  as  Echo 

was  from  Its  launch  on  12  August  1960,  to  24  August  1960,  the  perturbing 

force  due  to  solar  radiation  is  continuous  and  the  eclipse  factor  isV»  1 

in  the  expression  for  F  ;  The  integration  of  the  perturbation  equations 

0 

presents  no  problem  in  this  case.  More  often,  however,  the  satellite 
passes  through  the  Earth's  shadow  (where  V  *  0)  on  each  revolution.  The 
eclipse  thereby  produces  discontinuities  in  the  perturbing  force.  The 
customary  way  of  handling  this  discontinuity  is  to  evaluate  the  perturbative 
variations  in  the  elements  after  each  revolution.  The  values  of  the  satel¬ 
lite's  anomaly  at  the  points  where  the  satellite  leaves  and  enters  the  shadow 
on  each  revolution  are  then  the  lower  and  upper  limits  of  integration, 
respectively,  in  the  perturbation  equations.  The  quantities!},  to  ,  l  e,  and 
the  eclipse  points  are  considered  constant  during  the  one  revolution  interval. 
It  is  desirable  to  avoid  this  once-per-revolution  integration,  however,  and 
to  integrate  over  a  much  longer  time  interval. 

The  time  spent  in  eclipse  changes  continually  due  to  the  apparent 
motion  of  the  Sun  and  the  perturbations  on  the  satellite  orbit.  The  change 
per  revolution  is  most  rapid  when  the  satellite  orbit  plane  is  entering  or 
leaving  the  Earth's  shadow.  Once  in  the  shadow,  the  time  spent  iti  eclipse 
per  revolution  does  not  change  appreciably  over  a  number  of  days.  This  is 
illustrated  in  Figure  1,  which  shows  the  variation  of  time  in  sunlight  per 
revolution  for  Echo  around  24  August  1960,  when  it  was  ending  a  period  of 
some  12  days  of  continuous  exposure  to  sunlight.  The  change  was  10  percent 
the  first  day,  not  quite  5  percent  the  second  day,  and  it  quickly  decreased 
to  1  percent  per  day  and  less.  This  suggests  the  following  approximation, 
which  was  used  in  integrating  the  perturbation  equations  over  many  revolutions. 

The  parameters  were  integrated  over  many  revolutions  on  the  basis 
of  continuous  exposure  to  sunlight.  The  eclipse  factor  V  appearing  in  the 
expression  for  F^  was  then  given  the  form 


where  A  t  is  the  time  per  revolution  spent  in  the  shadow  where  the  perturbing 
force  does  not  act,  and  P  is  the  period.  During  the  intervals  when  Afc  is 
changing  very  slowly,  v  need  be  re-evaluated  only  every  day  or  so. 
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The  eclipse  factor  V  is  determined  according  to  the  following 
steps.  If  the  time  t  is  within  24  hours  of  the  time  t^  for  which  the  last 
eclipse  factor  V ^  was  computed,  this  section  is  bypassed. 

(a)  Calculate  the  geocentric  angular  distance  d  between  the  orbit 
plane  and  the  shadow  axis  (see  Figure  2): 


where 


W 


sin  (-  W  •  L  )  ■  sin 
—  — e 

Wx  ■  sin  sin  i 
W  *  -  cos  sin  i 

y 

W  ■  cos  i 


-1 


(-WL  -  W  L  -  W  L  ) 
x  xo  y  y®  z  z« 


xo 


yo 


zo 


cos  i 
e 

cos  €  sin  l 

s 

sin  €  sin  i 


The  angle  d  is  between  +  90  .  It  is  taken  as  positive  when  the  shadow 
axis  is  north  of  the  orbit  plane  and  negative  if  south  of  the  orbit  plane. 


(b) 


Calculate  the  angle 
(  7T-  -ty  )  -  sin  — 


If  |  7T  -  i/j  |  -  |  dg  |  <  0.01,  assume  V«  1,  that  is,  the  satellite  will 

not  pass  through  the  shadow,  and  omit  the  following  sequence  of  equations 
for  determining  V. 


If  I  tt  -  $  |  -  |  dj> 

the  satellite  enters, 
calculated. 


0.01,  the  eccentric  anomalies  at  the  points  where 

E,  ,  and  leaves,  E  ,  the  Earth's  shadow  must  be 
in  out 


-13- 


FIGURE  2.  ECLIPSE  GEOMETRY  PROJECTED  ONTO  CELESTIAL  SPHERE 


[ 
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(c)  The  values  of  the  eccentric  anomaly  at  the  eclipse  points  are  the 
solutions  of 

1  -  a2  (1  -  e  cos  E)2  +  a2  [(cos  E  e)  (P  »  ) 

+  *1  -  e2  sin  E(<}  •  L^)]  2  •  0  . 

There  are  only  two  valid  solutions  to  this  equation  since  the  bracketed 
quantity  must  be  negative  in  order  for  an  eclipse  to  take  place. 

Calculate  initial  estimates  for  E  at  the  beginning  (E^  )  and  end  (&out) 
of  eclipse: 

Ei„  ■  *  -  “ 

and 

A  u 

E  *  P  +  0  -<0 

out  2 

where 

A  u  ■  2  cos  *  j  — co—&  1  ;  A  u  <  180° 

!  cos  d  j  v 

l  a  J 

_L  C  •  M 

P  *  tan  ^  ^  ;  the  quadrant  is  determined  from  the  signs 

of  the  numerator  and  denominator. 

The  nodal  unit  vectors  N  and  M  are  shown  in  Figure  3.  Their  components 
are 

n  ■  cos  n 

x 

N  ■  sin  A 

y 

n  -  o 

z 

and 

M  ■  -  sin  Cl  cos  i 
x 

My  »  cos  (1  cos  i 
M  -  sin  i 


I 

I 
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FIGURE  3.  ORIENTATION  UNIT  VECTORS  AND  ANGLES 
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and  the  auxiliary  unit  vector  C  is  (cf.  Figure  2) 
C  »  -  sec  d  L  -  tan  d  W  . 


(d)  Calculate  the  dot  products  P  •  L  and  £  •  L  : 


*•*.•*{ 


(1  -  cos  i)  (1  +  cos  e  )  cos  (  ft  -(o  -  £  ) 

® 

(1  -  cos  i)  (1  -  cos  e  )  cos  (  ft  -  a>+  £  ) 


+  <1  +  cos  i)  (l  +  cos  e  )  COS  (ft  +  CD  ■  £  ) 

© 

+  (1  +  cos  i)  (1  -  cos  e  )  cos  ( ft  +  a>  -  £  ) 


•i 


+  2  sin  i  sin£  [cos  (co  -  £ m)  -  cos  (cd+  £ 

2  *  *  4  ■((!  "  cos  i)  (1  +  cos  e  )  sin  (ft  -  co  -  £g) 

+  (1  -  cos  i)  (1  -  cose  )  sin  (ft  -  0)  +  £) 

© 

(1  +  cos  i)  (1  +  cos  e  )  sin  (ft+  a>  -  £  ) 

-  (1+  cos  i)  (1  -  cos  e  )  sin  (ft+a)  +  i) 

I 

+  2  sin  i  sine  [-  sin  (  a)  -/0)  +  sin  +  [ 

(e)  With  the  initial  estimates  of  E^  and  E  ,  f(E)  and  f'(E)  are 
calculated  for  use  with  Newton's  approximation  met Sod.  Calculate 

f(E) 


1  -  a2  (1  -  e  cos  E)2  +  a2  [  (cos  E  -  e)  (P  *  L  ) 


and 


+  -  eQ2  sin  E  (fi  *  L^)  ] 


f*(E)  *  -  2  a  e  sin  E  (1  -  e  cos  E) 

2  J 

+  2  a  [  (cos  E-e)  (P’L)+V1- 


eQ  sin  E  (2  •  L#)] 


x  [ -  sin  E  (P  'L^)  +  ^  1  -  e2  cos  E  (g  *  L  )  ] 
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(£)  Calculate  the  correction  to  the  initial  estimates  of  E.  and  E 

in  out 


AE 


-  _L.CE) 

f'  (E) 


This  correction  is  applied  to  the  preceding  value,  E  1  »  E  +  AE. 
The  iteration  is  continued  until  *  ’ 


A  E 
E 


<  0.01 


(g)  The  eccentric  anomaly  at  the  eclipse  points  are  used  to  calculate 
the  mean  anomaly  at  these  points  from 


M.  -  E. 
in  in 


-  e  sin  E 


in 


and 


M  *  E  -  e  sin  E 
out  out  out 


These  in  turn  are  used  to  find  the  time  spent  in  the  shadow: 


A  t 


M 


out 


n 


from  which  the  eclipse  factor  V  is 


where 


V 


a 

o 


13  [  ! 

n 


|  sin2  io)  (l-K2)] 


2.2  EPHEMERIS  CALCULATION 


The  ephemeris  calculation  determines  a  satellite's  geocentric 

position  r  and  velocity  rata  given  time,  t,  from  the  initial  orbit 

elements  a  „  ,  a  „  ,  h  ,  L  ,  and  the  drag  coefficients  c  ,  and  d  ,  ut  an 
xN  yN  — o  o  o  o 

o  o 
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epoch  time  tQ.  The  perturbations  caused  by  the  Earth's  atmosphere,  by 
direct  solar  radiation  (when  significant),  and  by  the  asphericity  of  the 
Earth  (including  terms  due  to  the  second,  third,  and  fourth  zonal 
harmonics)  are  included  in  this  calculation. 


2.2.1  Preliminary  Computations 

The  following  quantities,  required  in  later  calculations,  are 
calculated  from  the  given  initial  conditions  (see  section  3.1.1). 

(a)  Compute  the  semi-latus  rectum  p  ,  the  eccentricity,  e  ,  the  mean 
semi-major  axis  aQ,  and  the  distance  to  perigee  q  : 


p  *  h  •  h 
o  — o  — o 


2  2  2 

h  +:•  h  +:-  h 
x  y  z 

o  o  o 


2  2  2 

e  *  a  +  a  „ 
o  xN  yN 

o  o 


a  -  p  /  (1  -  e  ) 
o  ro  o 

q  ■  a  (1  -  e  ) 
o  o  o 


(b)  Compute  the  three  orientation  angles:  iQ,  the  inclination; 

the  longitude  of  the  ascending  node;  and  Ot^,  the  argument  of  perigee: 


• I  * 

.  -1  o 

i  »  cos  /= — ; 

o  -/p; 


0  <  i  <TT 
—  o 


Cl 


tan 


m  tai, 

o 


-h 


yN, 


xN 


0 

of 

f 

o 

if 

f 

0 

if 

the  quadrant  is  determined  from  the  signs 

e  numerator  and  denominator,  0  <  Cl  < 

o 
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(c)  Compute  the  mean  argument  of  latitude,  U  ,  and  the  mean  value  of 

the  mean  motion,  n  :  ° 

o 

If  h  0  (direct  motion) 
o 


U 

o 


If  h  <0  (retrograde  motion) 
z 

o 


U 

o 


L 

o 


+  Q 

o 


n 


o 


(1 


where 


~  J,  a  2  -  1.623675  x  10'3  (earth  radii)2 

z  z  e 


and 


k  -  0.074  365  74 
e 


(earth  radii) 3^2  /  min 


(d)  Compute  the  drag  coefficients  c"  and  d: 


2 

360n  c 
o  o 


7T 

If  the  drag  acceleration  is  to  be  used,  d  is  read  from  storage  or,  if  zero 
is  found  there,  it  is  computed  from 


d  -  A  (c")2 


1  + 


3<nD  '  no) 
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where 


nD  =*  0.072  220  521  and  the 


empirical  coefficient  A  ■  8  (at  present)  ; 
if  (nD  -  nQ)  <  10  10,  set  (nD  -  nQ)  -  10  10 


2.2.2  Atmospheric  Dras 


Compute  the  purturbations  in  a,  e,  and  p  caused  by  the  Earth's 
atmosphere  during  the  time  interval  (t  -  tQ) .  At  time  t: 


a  *  a  1  +  2  c"  (t  -  t  )  +  3d  (t  •  t  ) 
o  o  o 


e  *  1  -  q  /a  for  a  >  q  and 
o  —  o 


e  =  0 


for  a  <  q 


p  »  a(l  -  e  ) 


2.2.3  Asphericity  of  the  Earth 

Using  the  initial  values  of  p  ,  n  ,  and  i  ,  calculate  the 
secular  bulge  effects  on  (2,  Cl)  ,  a^  and  a^. 

-i-£«fl+n+n  +  a  +  n 

dt  1  2  3  4  5 

a  2 

Q  3  T  e  ,  .x 

1  "  '  2  J2  “I  (C°S  io)  no 

Po 


3  T  2  e  .  9 

7T  J„  — y—  n  cos  l  7 

22  p  4  o  o[.4 

o 


3  v  l  -  e 


-  sin2  i  /5  9  /.  2| 

°  (2  2^  8o  >  J 
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*  -  77  J.  —^7  n  cos  i  (4-7  sin2  i  ) 
3  16  4  4  o  o  o 

po 


n  3  r  2  3e  2  J/1J.5.2.V 

4  8  J2  “ 4  eo  no  C0S  K  (1+  4  8ln  V 

po 


fy-  *  77J,  —^7  e  n  cos  1  (4-7  sin2  i  ) 

5  32  4  4oo  o  o 

Po 


n  ■  n  +  7-^  (t  -  t  > 
so  o  dt  o 


d  a) 

dT  "  “1  +  0)2  +  co3  +  +cd5 


“l  -  I  J2  V  (4  -  5  8in2  V  no 
Po 


/in  9  ...  2  e  , ,  ,-.2,v 

^2  "  4  J2  — T  (4  '  5  sin  V  no 

Po 


1  17  2  ,  2  ,  ,43  3  r.  2. 

1-  vl-e  -  sin  i  (75  -  fvl-e  ) 
O  O  4o  £.  O 


%  «  -  77  J,  -V  n  (16-62  sin2  i  +  49  sin4  i  ) 

3  32  4  4  o  o  o 

Po 

4 

“4  ■  if  J22  -V  %2 [  <4-5  *‘”2  V  <1+  s  ■1-2  V  -  f  -4  ‘„i 

Po 


-  -  775  J  -^7-  e  2  n  (24-84  sin2  i  +  63  sin4  i  ) 
D  LZ  o  4  4  o  o  O  O 
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W 


O)  * 

d  05 

(t 

t  ) 

s 

dt 

o' 

05  ■  05 

+ 

05 

SO 

0 

S 

w  - 

sin 

n 

sin  i 

X 

SO 

i 

w  - 

-  cos 

a 

sin 

cos  i 


z 

p 

l  „  * 

e 

(a 

COS  05  - 

a  » 

xNg 

e 

0 

^0 

S 

yN 

e 

(a 

sin  05  + 

a 

lyNs 

e 

0 

v  xN 

0 

s 

y*. 

sin  oo  ) 
s 

cos  05  ) 
s 


If  the  quantity  <  1.0  cm^/gm. ,  where  m  is  the  satellite's 

mass,  the  effect  of  direct  solar  radiation  on  the  satellite  orbit  is  not 
considered  and  the  ephemeris  calculation  proceeds  directly  from  Section  2.2.3 
to  Section  2.2.5. 


Ay  2 

If  the  quantity  ^  1.0  cm  /gm. ,  the  perturbative  effect  of 

direct  solar  radiation  on  the  satellite  orbit  is  computed  according  to  the 
formulation  outlined  in  Section  2.2.4.  This  limit  may  be  changed,  if 
experimentation  and  usage  show  that  the  radiation  should  be  applied  to  more 
or  fewer  satellites. 


2.2.4  Computation  of  Radiation  Pressure  Corrections 

The  computation  of  the  perturbative  variations  of  the  elements 
due  to  direct  solar  radiation  pressure  proceeds  essentially  as  in  Section  2.1. 
The  eclipse  factor,  v  ,  is  computed  as  previously  described,  except  for  the 
dot  products  (P  ■  L  )  and  (<J  •  L  )  which  are  reformulated  in  paragraphs  f 
and  g  below.  The  equations  for  IRe  corrections  to  the  parameters  e  cos  7T, 
e  s in  7T ,  Wx  >  Wy  >  and  L  have  also  been  rewritten,  as  shown  below,  for  pro¬ 
gramming  efficiency.  The  a  ^  coefficients  are  calculated  as  shown  in  Table  I 
using  the  relations  it  ■  _i_  and  05N  =■  1  dco  .  The  b  .  terms  are 

k  dt  k  dt 

6  6 

computed  as  the  reciprocals  of  the  a,  coefficients  with  il  ,  05  ,  and  l 

^  r  i  so  SO  9 

replacing  it,  C15N  ,  and  n@  respectively. 
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The  perturbations  in  the  elements  are  first  computed  at  t, 
eliminating  terms  whose  a.  coefficients  are  greater  than  10  .  The 
perturbations  are  then  recomputed  at  t^,  with  the  same  terms  eliminated. 

The  final  perturbations  in  the  elements  are  computed  as  follows 

A  (e  cos  7 r)  ■  A  (e  cos  id  -  A  (e  cos  7r)  tQ 

A  (e  sin  7 r)  ■  A  (e  sin  n)  -  A(  e  sin  7r)fco 


4WX  ■  (4Wx)t  - 

<4Wx>to 

4Wy  -  <4»y)t  - 

<aVto 

A  L  *  A  L  -  A  L 

t  to 


In  this  way,  the  very  long  period  terms  are  ignored;  the  moderately  long 
period  terms  are  treated  as  secular;  and  singularities  caused  by  very  long 
period  terms  are  eliminated. 

(a)  Computation  of  A(  e  CDs  7 r)  (DECPI) 

A  (e  cos  77)  -  kQ  |  kj  [  2k2  (-an  cos  bu  +  a12  COS  b12^ 

'  k3  (k4  a13  COS  b13  +  k5  S14  cos  b14}  +  k6  coa  lm 
+  k?  [  kg  (-au  cos  bn  +•  a12  cos  b^  +  «16  cos  b16 

'a17  C0S  b17) 

♦  k9  <“al4  cos  b14  "a18  COS  b18  +  a20  COS  b20} 

+  k1Q  <-a13  cos  b13  -  a19  cos  b^  +  a^  cos  b^) 

+  knl  cos  l  ] 
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(b)  Computation  of  A  (e  sin7T)  (DESPI) 

A  (e  sin7r)  =  kQ  |  [  2k2  (-a sin  b^  +  a12  sin  b12) 

^3  ^4  ai2  s^n  ^13  k5  ai4  8^n  bi4^  ^  k 26  8^n  .1 

+  k7  [k8  (_all  Sin  bll  +  3 12  Sln  b12  +  a16  Sin  b16. 

-  a17  Sin  b17> 

-I-  kQ  (-a^  sin  b,4  -aia  sin  b1R  -!  a,n  sin  b?f)) 


14  18 


:+k10  (-a13  sin  b13  -a19  sin  b19  +  a21  sin  bn> 

+  k30  sin  l%  ]  | 

(c)  Computation  of  Aw  (DELWX) 

AWx  -  kQ  k^2  [k^3  (-a^  sin  b^  +  a2  sin  b2> 

-!  k14  (a3  sin  b3  *  a4  sin  b4> 


k15  a5 

sin  b^  -i 

k16 

a6 

sin  b^  -!  k^y 

a7 

sin  by 

+  k18  a8 

sin  bg  - 

kl9 

a9 

sin  b9  -k2Q 

aio 

sin  b1Q  ] 

Computation  of  AWy 

(DELWY) 

AWy-koki2  [k13  (a 

^  cos  b^ 

-a2 

COS 

V 

+  kl4  (a3 

cos  b3 

'+  a4 

cos 

V 

+  k28  a5 

cos  b^ 

k29 

a6 

cos  b.  -kn_ 
o  1/ 

a^  cos  by 

*18  a8 

cos  b  -! 

O 

k19 

a9 

cos  b9  k20 

aio 

cos  bio  ] 
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(e)  Computation  of  A  L 
L 


(DLTL) 


ko  {k21  C  k22  al  8in  bl  +  k23  a2  Sin  b2  +  k24  a3  8in  b3 
+  k25  a4  sin  b4  +  k2  (-a5  sin  b$  +  afe  sin  b&)  ] 

+  k27  [k1Q  (-aL  sin  bx  +  a3  sin  b3  ) 


+  kg  (a4  sin  b4  -a2  sin  b2> 

+  kQ  (-a  sin  b  +  a,  sin  b,)  ] 


(f)  Computation  of  (P  •  L  )  (PD0TL) 

—  —8 


4(P  •  LJ  =  k25cos 


'K1 


+  k24e0*(b3l+  k23co,lb2U 


+  k„„cos  jbj^j  +  k2  cos  jb^^j  -  cos  jb^JJ 


22 


(g)  Computation  of  (Q  •  L  )  (QD0TL) 

—e 


4((J  '  L-)  *  k2s  sin  +  k_4  sin 


-  k22  sin  b1  +  k2 


31]  '  k23  8in  (hi') 


-  sin  |  b^1  +  sin/b  * 
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l 

mr 

TABLE  II 

i 

The  Coefficients 

k  -  k 

K0  30 

1“  Coefficient 

Symbol 

Formula 

L 

XYZKO 

= 

-  —  [PSUN-AG0M-XYZNU] 

[ 

XYZK1 

= 

RTAUP • RTE0SQ 

r  k2 

XYZK2 

= 

2- SINI- SINEP 

L  k3 

XYZK3 

= 

2(1  -  C0SI) 

r  k4 

XYZK4 

= 

(1  -  C0SEP) 

1  k 

R5 

XYZK5 

= 

(1  +  C0SEP) 

[ 

XYZK6 

= 

XNSDN  [<1  +  C0SI)  (C"SEP)1 

k7 

r 

XYZK7 

= 

RTAUP  (Eli)2  SINI 

RTE0SQ  (1  +  C0SI) 

l  k8 

XYZK8 

= 

2  C0SI- SINEP 

r  k9 

XYZK9 

= 

SINI  (1  +  C0SEP) 

t-  k 

10 

XYZK10 

= 

SINI  (1  -  C0SEP) 

[  kll 

XYZK11 

= 

XNSUN  (SINI)  «*“» 

k12 

XYZK12 

= 

RTAUP 

RTE0SQ  E  ^ 

L  k 

13 

XYZK13 

= 

2  C0SI  SINEP  (1  +  C0SI) 

[  k14 

XYZK14 

= 

2  C0SI  SINEP  (1  -  C0SI) 

i— . 

k15 

XYZK15 

= 

2  SINI  (C0SI  +  C0SEP) 

[  k16 

[ 

XYZK16 

2  SINI  (C0SI  -  C0SEP) 

l 

[ 

[ 

-27- 

I 

I 

I 

I 
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[ 
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[ 

[ 
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[ 


[ 

l 

l 

I 

[ 
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(h)  Application  of  Perturbations  In  the  Elements 

The  corrective  quantities  A  a,  A  (e  cos  7r),A(e  sin7r),  AW  » 

AW  ,  and  AL  due  to  direct  solar  radiation  are  now  applied  to  the  orbital 
elements  at  time  t.  The  perturbative  effects  of  the  Earth's  asphericity 
and  atmosphere  have  previously  been  included  in  the  elements  a,  a^  »  e  cos  4), 
*  e  sin  oo,  ft  ,  and  L.  Since  the  radiation  pressure  perturbations  are 
based  instead  on  the  parameters  (e  cos  tt),  (e  sin  7r),  W  ,  and  W  ,  they  must 
be  introduced  in  the  following  manner.  X  ^ 

(1)  Calculate  (e  cos  tt)  and  (  e  sin7T)  at  time  t: 

e  cos  i  •  [a  cos  !1-  a  sin  !)]+  A(e  cos  tt) 

xN  yN 

e  sin  tt  *  [a  „  sin  ft  +  a  „  cos  ft  ]  +  A  (e  sin  tt) 
xN  yN 

The  square-bracketed  terms  include  the  bulge  and  drag  perturbations  at  time  t. 


(2)  The  eccentricity,  which  now  includes  the  perturbative  effects 
of  the  bulge,  atmospheric  drag,  and  the  radiation  pressure,  is  computed  from 

e  *  [  (e  cos  tt  )2+  (e  sin  tt)2  ]  1^2 

(3)  Similarly,  W  is  computed  from 
W  *  [W  ]  +  AW 

X  X  X 


W  -  [  W  ]  +  A  W 

y  y  y 


[»,)  - 


W  AW  +  W  AW 

X  X  V  v 

cos  i 


where  the  square-bracketed  components  of  W,  computed  in  2.2.3  above,  include 
the  bulge  perturbations.  If  j  cos  i  |  is  too  small,  compute 

2  o 

|  W  I  *  /  1  -  W  -  W  and  apply  the  sign  of  [-  W  AW  -  W  AW 

z  J  x  y  rr/  °  x  x  y  y 
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(4)  The  inclination  i  and  ft  are  now  revised,  both  radiation  pres¬ 
sure  and  secular  bulggi  effects  being  included: 

sin  i  ■  ^  1 


cos  i  ■  W 


- j— 

cos  i 


i  ■  tan 
-  W 


-1 


sin  i 


cos  i 


cos 


SI 


;  0  <  i  <  180c 
W 


-r^-r  ,  sin  A  ■  —~~7 

sin  i  sin  i 


A  »  tan 


-1 


sinfi 
cos  A 


;  0  <  A  <  2  7T 


(5)  Recompute  a^,  a^N  with  the  new  SI  • 

a^  ■  (e  cos  7T )  cos  (ft  +  (e  sin  7t)  sin  SI 


yN 


(e  cos  7 r)  sin  ft-f-  (e  sin  7T)  cos  ft: 


(e  cos  7r)  and  (e  sin  7T)  are  given  by  2.2.4  (h),  (1). 


2.2.5  Geocentric  Position  and  Velocity  Calculations 

(a)  Calculate  the  long-period  effects  on  ft  ,  a^,  and  ayN’ 

ft  *  A,+'  A_4r  ftQ  +  ft.  +  ftin 
L  6  7  o  9  10 

1  ^  7  (7_15  811,2  U 

*  A  =  -  -  J  — =—  e  cos  i  - x -  sin  2  0) 

6  8  2  o  c  .  2  ,  v  so 

p  (4-5  sin  io) 

a  2  (14-15  sin2  i  ) 

**ft7  -  -  —  J  —  e  sin  i  cos  i  - 5 - r  sin  2  0) 

7  16  2  2  o  o  ,,  ,  .  2  .2  so 

p  (4-5  sin  i  ) 


A, 


.  J_  a  cos  i 

1  3  e  o 

-  —  —  e  — — 7-  cos  (£> 

2  J 2  P  sin  i  so 


Note:  All  asterisks  denote  program  size  tests  explained  at  the  end  of  this 
Section,  2.2.5 
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where 


1  J3  -3 

i  ~r  a  =  -  1.1548  X  10 
z  ^2  e 


.  J.  a  2  (3-7  sin2  1  ) 

^  o  5  4  e  2  .  o 

*  V a  -  '  a  ~r  ~o  e  cos  i  - r -  sin  2  (X> 

9  oj-2  o//c.2.v  so 

2  p  (4-5  sin  i  ) 


where 

4  —  a  2  =  -  1.0682  x  10‘3 

8  J2  e 

J,  a  _  (6-7  sin2  i  ) 

**  "  -  —  — - re  cos  i  sin  i  - r — —  sin  2  CD 

10  16  J2  p2  °  °  (4-5  sin2  io) 

fl  >  (2  +  £2 

sL  so  L 

a  *  a  "fa  +  a  +  a  -ha  +  a  +  a  -fa.  -fa  -fa  -fa 
xNl  xN][  xN2  xN3  xN4  xN5  xN6  xN?  xNg  xN9  x*l0 


l  ae2  2  (14-15  sin2  i  ) 

*  aw  *77  T  esin  i  r 

1  p  (4-5  sin  i  ) 


cos  0) 


so 


,  a  2  Q  (13-15  sin2  i  -fr  105/2  sin^  i  ) 
lei  o  o 

*  a  m  -  —  — —  e  . -  cos  CD 

xN  32  2  ..  2  .  so 

2  p  (4-5  sin  i  ) 


Note:  All  asterisks  denote  program  size  tests  explained  at  the  end  of  this 
Section,  2.2.5 
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,  a  2  ,  ,  (13-15  sin2  i  )  (14-15  ain2  1  ) 

**  J,  '  --f  e3  sin2  1  - J - r - S-  cosOJ 

xN.  32  2  2  O  //  c  j  ^  i  \" 

3  p  (4-5  8 in  iQ) 


1  ae  3 

*  axN  -ah’?  6 

4  p 


(14-93  sin2  i  4'  sin4  i  ) 
' _ o  2 _ o_ 

(4-5  sin2  i  ) 
o 


cos  3co 


1  a  2  (13-15  sin2  i  )  (14-15  sin2  i  ) 

**  a  =—  J0  -f  e  sin  i  - - - 5 - 5 - 2  cos  3  CD 

xNc  32  2  2  o  //  c  .t  2  *  \2  80 

5  p  (4-5  sin  i  ) 


^  _  1  J3  %  2  2  (1-2  sin  ip) 

*xN  4  J-  p  6  S  n  o  sin  i  8  n  ^so 

6  2  o 


_  2  (6-7  sin  i  ) 

e  8  0  o' 

5  4  e  .  2  . 

*a„«—  —  — re  sin  i  - 5 -  cos  m 

**7  16  J2  P2  °  (4-5  sin2  i  )  60 

o 

c  J,  a  2  _  (6-29  sin2  i  ■+•  49/2  sin4  i  ) 

5  4  e  3  o  o 

a  . _  -  -  —  —  e  - ^ -  cos  a) 

**8  32  J2  p2  (4-5  sin2  i  )  80 

o 


5  J4ae23  2  (6-7  sin2  i) (13-15  sin2  i) 

**  3  W  =  ^  T  9  e  sin  i  - 2 - 5 - 5 -  COS  CD 

xNft  32  2  c  .  2  N2  so 

9  2  p  o  (4-5  sin  i  ) 

o 


o  o  77  4 

c  J.  8  ,  (6-41  sin  i  +  sin  i  ) 

5  4  e  J  o2  o 

*  *xb,„  ■111  1  - -  77 t  ;-. - c°* 3  so. 

10  2  p  (4-5  sin  ip) 


Note:  All  asterisks  denote  program  size  tests  explained  at  the  end  of  this 
Section,  2.2.5 
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s  J4  a  3  2  (6-7  sin2  i) (13-15  sin2  i  ) 

**  a  “  f2  “  -f  e3  sin2  i  - 2— 5 - 5 - 2_  cos  3  o> 

^11  32  J2  p2  0  (4-5  sin2  i  )2  80 

o 


a  XT  =  3  XT  ^  3  M  +  a  M  a  „  +  3  XT  +  a  XT  +  3  XT  +  3  *T  +  3  „ 

yNL  yNj  yN2  yN3  yN4  yN5  yNg  yN?  yNg  yNg 


+  a  :+  a  +  a  +  a 
yN^  yNu  yN12  yN13 


,  a  2  2  (14-15  sin2  i  ) 

*  a«  “  “  T7  Jo  ~f  esin  i  - 5 - —  sin  CD 

yNl  16  2  p2  0  (4-5  sin2  i  ) 

o 


2  «  ,  2  ,  .  105  ,  4  .  . 

,  a  (14-65  sin  i  +  — —  sin  1  ) 

1  „  e  3  o  2  o 

*a  «— -J  -^e  - - -  sin 

yN  32  2  2  2  .  so 

2  p  (4-5  sin  i  ) 

o 


,  a  2  .  ,  (13-15  Sin2  i  )  (14-15  sin2  i  ) 

•JUJU  --  T  ®  ^  f  ^  »  O  O 

*yN  "  32  J2  ~2  *  8l"  *o  - -  7“  2  7  ,2 -  sl”  “so 

7  3  p  (4-5  sin  1  ) 

o 

,  a  ^  „  (14-93  sin^  i  +  sin^  i  ) 

.  1  _  e  J  o  2  o  .  _ 

3  <jo  n  ®  9  S  m  3  (I) 

yN4  32  2  2  (4-5  sin2  i  )  so 

o 


,  a  ‘  (13-15  sin  i  ) (14-15  sin  i  ) 

**  a  »  —  J  -f  e  sin  i - 2 - r - 5 - —  sin  3gd 

yNc  32  2  2  o  c  .  2  .  .2  so 

J  5  p  (4-5  sin  1  ) 

o 


1  a 

13  e 

a  =  -  T  —  —  sin  i 

yNg  2  J2  p  o 


Note::  All  asterisks  denote  program  size  tests  explained -at  the  end  of  this 
Section,  2.2.5 
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.  J-  a  2 

1  _ 3  _ e  e 

yN?  4  J2  p  sin  iQ 


i  J~  a  -i  (1-2  sin  i  ) 

a  -  T  —  —  e2  - - — ~~  cos  2  CD 

yN„  4  J  p  sin  i  so 

O  ^  O 


c  j,  a  2  ,  (6-7  sin2  i  ) 

*  a  *  -  -r,  —  — f  e  sin  i  - - — —  sin  CD 

yN9  16  ;J2  2  o  (4_5  sin2  } 

o 


2  2  4Q  4 

J.  a  ~  (6-29  sin  i  +  —  sin  i  ) 

*  aN  -f  e3  - - o  sin  m 

yN10  32  J2  p2  (4-5  sin2  i  )  8 

o 


c  J.  a  2  ,  ,  (6-7  sin2  i  )  (13-15  sin2  i  ) 

**  a  M  ~  -J  e3  Sin2  i  - 2 - s - 5 - *  sin  CD 

yNll  32  J2  p2  °  (4-5  sin2  i  )2 

o 


2  2  77  4 

,  J,  a  (6-41  sin  i  +:-  —  sin  i  ) 

5  4  e  3  o  2  o  . 

«  —  ■—*  — —  e  - * -  si 


ayN19  ■  32  J  2  e 


sin  3  CD 


(4-5  sin  i  ) 
o 


J,  a  2  (6-7  sin2  i) (13-15  sin2  i  ) 

......  54e3.Z.  o  o  .  „ 

**  a  „  =777—  -7  e  sin  1  - 5 - — - sin  3  0) 

yN. -  32  J,  2  o  ,,  c  .  z  2  so 

J  13  2  p  (4-5  sin  i  ) 

r  o 


a  ■  a  4.  a 

^  SL 


Note:  All  asterisks  denote  program  size  tests  explained  at  the  end  of  this 
Section,  2.2.5 
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yN 


SL 


a  +  a 
yNS  yNI 


(b)  Compute  the  mean  longitude  L  at  time  t: 

L  -  L  +  n  (1+AM+A7T)  (t-t  )  +  n  c"  (t-t  )2  +  n  d  (t-t  )3  +  L.  +  AL 

OO  O  O  o  o  O  li 

where  n  AM,  n  Att  ,  and  L  are  computed  as  in  the  following  three  subsections, 
and  AL  fs  the  effect  of  solar  radiation  pressure  (Section  2.2.4)  or  is  set 
equal  to  zero  if  radiation  pressure  effects  are  not  considered. 

(1)  First  and  second-order  secular  changes  in  the  mean  anomaly  M: 

n  AH  -  M,  +  M.  +  M, 
o  12  3 


M  3  _  2  ae  2  25  5  nr- 

M1  ■  2  J2  4  J  1_eo  32  2  J  l* 


+  (||  -  3  i-e  2)  sin2  i  +  ( 
lo  ^  o  o 


9  n  4  1 

32  "  2  J  1_eo  )C0S  ioJ  "o 


m  15  I  2  I. 
M2  "  64  J2  4  v/1” 
P, 


2  .13  .  .  2  .  5  4  . 

a  (“T  -  9  sin  i  -  —  cos  i  )  n 
o  2  o  2  o  o 


45  e 

M3-128J4  T 

Po 


1-e  2  e  2  (27-30  sin2  i  -  35  cos^  i  )  n 
o  o  o  o  o 


(2)  First  and  second-order  secular  changes  in  the  argument  of 
perigee  7T  : 


n  A7T 
o 


W1  +  \  +  w3  +  \  +  ^ 
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TT  .  A  J 


12  2  2 

p„ 


e  (2  -  ~r  sin 2  i  -  |  cos  1  |  )  n 


o  o 


7T  .  2  j  2 
2  4  J2 


\|(4-5  sl„2l0)  [  l-  7  i-«o2  -  (£§  -  I  «4T2).m2lo  ] 
p0  L 


-  I  COS  iD|  Cl  '  2 


3  “  J  1-e 


2  -  sin2i  (|  -  3  v  1-e  2)  ] 
o  o  '3  o  ^  I  o 


9  .  2  e  2 


J,  -7  «.  j  (4-5  sin  ij  d+^4  8in  i0)  -  %  cos^i 


7T  _  2 

3  16  “2  4  o 

Po 


1  „_2.  v  5  4. 


-  7  |  cos  i  |  (4+5  sin2i 


'In 

o  r  o 


\  ■  ’  M  J4  ~h  C  16'62  811,2  ^  +  49  811,4  lo 
Po 


-2  I  cos  i  (4-7  sin  i  )  ]  n 
i  o  o  4  o 


^5  *  '  128"  J4  eo  t  72-252  811,2  10  +  189  8ln\ 

Po 


-12  |  cos  i  |(4*7  sin  i  )  ]  n 
1  o'  o  o 
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(3)  Long-period  changes  in  L: 

LL  “  L1  +  L2  +  L3  +  L4+  L5  +  L6 


*  L  ae2  J 

1  16  2  p2<4-5  .l„2t„)  | 


[<l-e2)3/2  -l  J  (14-15  .lo2l  )  .lt,2l 

o  o 


+  ^  [14-79  sin2i0+  ~  sin4iQ  -  |  cos  i  |  (14-30  sin2i  )  ]  Isin  2 


,  a  sin  i  (14-15  sin  i  ) 

**  t  .  —L  j  l  _ o  _ o' 

2  16  J2  2  2 2 
P  (4-5  sin2  i  >Z 

o 


[  (13-15  sini 


5  I  cos  i  I  ]  sin  2  a) 
o  1  so 


Il.i  i  !e  [ 

'  3  2  J„  p  L  e 


sin  i  ]  cos  U) 
o  so 


1  3  e  |  cos  i  | 

L4mi  t2  r  e  Tirf  ( i cos  i0i  _1>  cos  “s 


*  L5“  “  h  p2(4-5esin2i  )  ' 

o 


I  [(1-e2) 


3^-l]  (6-7  sin2i  )  sin2i 

o  o 


+  e“  [6-35  sin2 i  +  sin4i  -  i  cos  i  I  (6-14  sin2i  ) ]  Isin  2  0) 

1  o  O'  o  J  f  so 


Note:  All  asterisks  and  daggers  denote  program  size  tests  explained  at  the 
end  of  this  Section,  2.2.5 
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** 


5  J4 

2 

a 

e 

2 

sln^i  (6-7  sin^i  ) 
o  o 

16  J2 

2 

P 

e 

(4-5  sin^i  )^ 

0 

-5 

]  cos 

i  1 

]  s  in  2  CD 

0 

so 

Compute 

the  mean 

argument  of  latitude. 

U  »  L  - 

flsL 

if 

W  >  0 

2 

U  =*  L  + 

flsL 

if 

W  <  0 

2 

[ (13- 15  sin  iQ) 


(4),  Solve  the  following  modified  form  of  Kepler's  equation  for  the 
quantity  (Eh CD)  by  iteration,  using  U  (mod  £n)  as  a  first  guess 

E+cu  -  U  +  a„  sin  (Ed-CD)  -  a  cos  (E+cd) 

^SL  yNSL 

(e)  Compute  the  geocentric  position,  r,  and  velocity,  r,  at  time  t 
by  means  of  the  following  sequence  of  equations: 

e  cos  E  ■  a  N  cos  (E+CD  )  +  a  „  sin  (E+cd  ) 

SL  7  SL 

e  sin  E  «  a  „  sin  (E-HD  )  -  a  „  cos  (E-+- CD) 

*”SL  7  SL 

2  2,2 

eT  *  a  „  +  a 

L  ^SL  yNSL 


r  *  a  (1-e  cos  E) 


s/'ua 

r  ■  -C1— 


e  sin  E;  ■  1 


Note:  All  asterisks  denote  program  size  tests  explained  at  the 
end  of  this  Section,  2.2.5 
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r2  “  4  J2  if  8in2ioL  COf  2  u 

li 


1t”r3  *  "  2  J2  ^  -  | 


(1 


OF-) 


(4)  Short- period  terms  in  r: 


72  ‘  2  'Al? 

Pi  L  “■  e_ 


3/2 


I  (i-7iT;7>]  < 


a  sin  u  -  a  cos  u) 
^SL  yMSL 


1tf2  "  2  J2  V  [¥ 


J2  J  £  U  '  *  *  Al> 

PL  L 


[ 


(±E) 


* 1  (Hl'  *4>,1°  2u  ' 2  Hl  Ht  “* 2u 


] 


Mote:  All  daggers  denote  program  size  testa  explained  at  the  end 
of  this  Section,  2.2.5 
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It  r3  “  +  8  J2  _  2 


1- J  1-e  "  ,  , 

- = — —  (1-4  sin  i  )  (a  sin  u-a  cos  u) 

PT  -  2  '2  oL  xNcl  yNs 


L  eT 


SL 


SL 


^3e2  -4 


x  1  3e  -4(a  sin  u-a  „  cos  u) 
**SL  y  SL 


!] 


r,  -  -  t  J-  sin  i  T  a  „  a 

4  4  2  pT  2  \/  PL  oL  ^SL  ^SL 

L 


r5  = 


-  '  T  J 


f1 

i  Pt 


2  -2  2  ,  pL  -Vfe-  ”> 


*6  =  '  \  J2  “*2  /  £  Sin  ^L81"  2U 
PL  / 


r7  = 


—  J 
8  2 


2  1 - 

!£_  flL 

2  '  Pt 
PL  L 


2  2 

eT  sin  i  T  sin  2u 
p  L  oL 


r8  = 


-  -  4  J. 


“%T  r?  sin2ioL^axN  Sin  3u_a«*>  cos  3u^ 


2  2  Pt  2 


SL 


yN 


SL 


r9  = 


8  J2 


2 

ae  ru  .2.  r  2  o  2. 

—  /  *"  loL  (a*»„  -  > 

PL  1  PL  L  SL  SL 


sin  4u-2a  „  a  „  cos  4u 
^SL  ^SL 


(5)  Short-period  terms  in  rv  (r  v 


r2v  »  J  up): 


(rv)  =  rv^  +  rv2  +  rv3  +  rv^  +  rv^  +  rvg  +  rv7  +  rvg 


•  rv  a 

rv^  *  -  —  Ar 


Note:  All  daggers  denote  program  size  tests  expiated  at  the  end  of  this 
Section,  2,2.5 
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rv 


2  4 


!  j2  \  [J  F\  <i-!*i»2w  <t> 


»L  '  'L 


e  /  M  _i_2 


™3  "  4  J2  2  V  P 
PL  L 


rv.  *  zr  Jn 
4  2  2 


!*!  I  JL 
PL2  J  h. 


sin  1  _ cos  2u 
oL 


2  2 

eT  sin  i  _  cos  2u 
L  oL 


rv5«-  J2 


ae  Tm. 

PT2  iPL 


sin  i  (a  cos  u-a  sin  u) 
oL  ^SL  ^SL 


.  5  T  ae  f  Jl_ 

rv6  =  8  J2  ~2  £ 
PL  si 


sin  i  T (a  „  cos  3u+a  sin  3u) 
°L  ^SL  yNSL 


**7  "  8  J2 


ge  /  M_ 

pl2  J  PL 


2  2  2 

sin  i  (a‘  -  a  ) 

oL  ^SL  ^SL 


™8  *  8  J2 


JL_  I  J£ 

PL2  '  PL 


2  r  2  2 

sin  i„T  I  (a  -a  )  cos  4u+2a  „  a  sin 

oL  L  ^SL  ^SL  ^SL^SL 


(6)  Short-period  terms  in  ft  : 


Aft  =  ft  u  +  ft12  +  ftu  +  ft14 
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3 

nil  =  4  J2~h  COS  1oL  sin  2u 
PL 


3  °g 

^12  =  -  2  J2  cos  \>L  (u-lHaxNgL8in  “'V^08  U> 


fi13  =  !  J2  -%  cos  ioL^3xN  8in  u  +  ayHCTCO“  u) 
p_  Sli  bL 

L 


1  de 

^ =7  -^T  COS  i  (a  sin  3u  -  a  cos  3u) 

14  4  2  2  oL  xKL.  yN 


(7)  Short-period  terms  in  i: 


A  i  =  i4  +  i5  +  i6 


H  =  f  J2  "^2  Sin  2ioL  C0S  2u 
PL 


i5  =  I  J2  Sin  2ioL  (axNQTCOS  U  -  ayN  Sin  u) 

p  OL  Oli 

Li 


h  =  8  J2  -h  Sin  2ioL  (axN  C0S  3U  +  ayN„Sin  3U) 

SL  dll 


SfbirdJ^gtcr^om/icuny-, 

AERONUTRONIC  DIVISION 


(8)  Orientation  vectors  W,  M,  N: 


i.  -  i  T  +  Ai 
k  oL 


cl  >  ci  +  a ci 

k  aL 


W  *  sin  f!,  sin  i, 
x  k  k 


W  =  -  cos  (2.  sin  i, 
y  k  k 


W  =  cos  i. 
z  k 


M  *  sin  ft.  cos  i, 
x  k  k 


M  =  cos  ft.  cos  i, 
y  k  k 


M  =  sin  i. 
z  k 


N  =  cos  ft  . 
x  k 


N  *  sin  ft , 

y  * 


N  -  0 
z 


/ 

\ 


(9)  Direction  vectors  U,  V: 


u.  =  u  +  A  u 
k 


U  =  N  cos  u.  +  M  sin  u, 
-  -  k  -  k 


V  =  -N  sin  +  M  cos 
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(10)  Position  and  velocity  r  and  t; 


rk  «  r  +  Ar 


i-k  -  r  +  A  r 


(rv)k  =  rv  +  A  (rv) 


r  *  r.  U 
—  k^- 


r  »  r^  +  (rv) k  V 


Size  Tests 


*  If  (4-5  sin2  iQ)  <  10  ^  disregard  tern 


**  If  (4-5  sin2  iQ) 2  < 10  ^  disregard  term 


r 


_3 

If  e  <10  expand  the  quantity 


It 


3/2 

(1-e2)  -1  3  .  1  2. 

e -  -  -  ^  e  (1  -  -  e  ) 


-3 

If  e  <  io  expand  the  quantity 


1  ’  ~  1'CL  =  1  1  2 

2  2  ^  +  4  el? 


t 


-3 

If  e  <  10  expand  the  quantity 
„  3/2 


1  -  U-V 


3  1  2. 

2  (1  "  4  elP 
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2.2.6  Deviations  from  the  Nominal 


Determine  the  radial,  transverse,  and  orthogonal  components  of 
the  deviation  in  position  from  mean  elements  for  the  nominal  case  (including 
all  terms) : 

radial  Ar 

transverse  r  A  9^ 

orthogonal  rA  9^ 

I  Ar| 

2.2.7  Subsatellite  Point 


=  rAu  +  r  cos  i  T  Aft 

Oh 

=  r  sin  u  A  i  -  r  sin  i  ,  cos  u  Aft 

oL 


=  [(Ar)2  +  r2  (A03)2  +  r2  (  AQJ2  ] 


1/2 


As  each  point  is  calculated,  the  position  obtained  is  used  in  the 
subsatellite  point  computation  for  output. 


*3  0  o  "1  O  / 

h=r-l+(-rf  +f)U  -  r  U  converted  to  kilometers 

im  Z  Z  Z 


u 


Latitude  =  tan 


/K  d-f)‘ 


converted  to  degrees 


Longitude  =  tan  ^  (y/x)  -  9^  -  9  (t- tQ)  converted  to  degrees 

•  -3 

where  9  =  4.375269511x10  radians  per  minute 


and  9  -9  +  0.9856472D  +  360.9856472F  converted  to  radians 


gro  gr00 


where  9 


is  the  Greenwich  sidereal  time  at  the  start  of  the  epoch 
oo  year  in  degrees,  and  D  and  F  are  the  days  and  fractions 
of  a  day,  respectively,  from  the  start  of  the  epoch  year 
to  the  epoch. 
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2.2.8  Position-Error  Analysis 

The  position-error  analysis  applies  the  above  formulas  2.2.6 
to  the  computed  differences  between  the  nominal  orbit  values  and  the 
experimental  (omitted  term) case.  The  subscript  n  refers  to  nominal  case 
quantities. 


Ar  =  r,  -  r, 
kn  k 


Au  = 


“kn  '  Uk 


a  =  ft,  -  n 

kn  k 


A  i  =  i 


^kn  Xk 


rA0=r  Au+r  cos  i  T  Aft 
3  n  n  oL 


r  A0,  =  r  sin  u  Ai  -  r  sin  i  _  cos  u  ao 
In  n  “  n  oL  n 


I'Ar  |  = 


(  Ar)2  +  r2  (A03)2  +  r2  (A 6 


•>*] 


1/2 


2.3  DIFFERENTIAL  CORRECTION  FORMULATION 

The  purpose  of  the  differential  correction  is  to  relate  the  topo- 
centric  observation  residuals  to  improvements  in  the  orbit  parameters.  The 
formulation  used  in  this  program  is  similar  (except  for  minor  modifications) 
to  that  appearing  in  Aeronutronic  Publication  U-880*,  and  is  presented  here 
for  reference. 

2.3.1  Compute  0^  ,  the  Epoch  Greenwich  Sidereal  Time: 

0  (deg)  =  (  0  -  360)  D  +  0F  +  0 
gro  8  oo 

*  Astrodynamic  Analysis  for  the  National  Space  Surveillance  Control  Center, 
Sections  5.1  and  5.2  of  Appendix  4A,  June  1,  1960. 
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Where  9  is  the  rotation  rate  of  the  Earth  in  deg/solar  day,  9 

*roo 

is  the  Greenwich  sidereal  time  at  the  start  of  the  epoch  year  in  degrees,  and 
D  and  F  are  the  days  and  fractions  of  a  day,  respectively,  from  the  start  of 
the  epoch  year  to  the  epoch. 


2.3.2  Station  Vector 


Compute  the  station  vector  R  (X.Y,Z,)  from  the  given  quantities, 
0,  the  geodetic  latitude;  the  east  longitude;  H,  the  altitude;  and  t, 
the  time  of  observation: 


C  = 


^  1  -  (e2sin20  ) 


2  2 
e  =  2f  -  r 


where  f  =  ^  is  the  Earth's  flattening 


S  =  C  (1-e  ) 


9  =  9  (t-t  )  +  9  +  X_. 

(rad)  o  gr  E 

o 


where  9  is  the  rotation  rate  of  the  Earth  in  radians/solar  minute. 


X  =  -  (C+H)  cos  0  cos  9 
Y  =  -  (C+H)  cos  0  sin  9 
Z  =  -  (S+H)  sin  0 

H  has  units  of  equatorial  Earth  radii 
where  X,  Y,  and  7  are  the  components  of  the  station  vector  R. 


2.3.3  Compute  the  Partial  Differential  Coefficients  R^  and  IK: 
2 

(a)  R  =  — —  esinE 
u  r 
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(b) 

H 

*C 

(c) 

R  - 
xN 

2 

8  r 
r  t  axN 

(d) 

R  = 
yN 

2 

a  r 
r  ^  ayN 

2R 

(e)  R  m  (U-U  ) 

c  n  o 

3Rn  2 

(f)  Rd-^  (U-Uo)2 

n 


(g) 


2  - 

U  =  —  Jl-t 
u  r  » 


(h)  U  =  (U-U  )  U 

n  o  u 

2  / 

(i)  UxN=  —•  \(1  +  ^  )  8in  (E  +  oj)  +  e  sin  E 


'-«/  1-e2  (1  +  J  1-e  )  2  .  1+  J l-eZ 


(j)  UyN=  —  ^  -  (1  +  ^  )  cos  (E  +a>)  +  e  sin  E 


■  [■ 


e2  -  a 


+  J 1-e2) 


1 


e  cos  E 


*/l-e2  (1  +  */  1-e2)  2 


1  +  'Z  l-e^  : 

J 
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E 

C 

C 

0 

c 

0 

L: 

[ 

[ 

i 

L 

l 

E 
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2U 


(k) 

u 

_ n 

(U  - 

c 

n 

3U 

(1) 

ud 

n 

2 

n 

<u  - 

2.3.4  Topocentric  Position  of  the  Satellite 

Compute  the  topocentric  position  of  the  satellite  at  the  given  time 
from  the  geocentric  position  £  and  the  station  position  vector  R: 


(a) 

fic 

(b) 

Pc 

(c) 

L 
— c 

=  ££ 
c  Pc 

where  the  subscript  c  denotes  computed  quantities. 

2.3.5  Slant  Range  Observations 

If  p,  the  slant  range,  is  observed,  compute  the  residual 
AP  =  Pq  -  P  ,  where  the  o  subscript  denotes  observed  quantities,  and  form 

the  partial  differential  coefficients: 


An  _  L  -  UR  +  L  •  V  U 

-  -  — c  —  n  — c  —  n 

n 


C.  =  L  •  UR„+L  •  V  U  „ 
ia  „  — c  —  xN  — c  —  xN 

xN 
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C  ,  -  L  '  U  R  „  +  L  •  V  U  „ 

Aayl|  -c  -  yH  -c  -  yH 


:A„  -  L  -  UR  +  L  .  VB 

A  U  — c  —  u  — c  —  u 

o 


C  a  r,  *  L  •  V  r  cos  i  -  L  •  W  r  sin  i  cos  u 
Aii  — c  —  — c  — 


C  a  ,  *  L  •  W  r  sin  u 

A  i  — c  — 


Ca„«L  -  UR+L  .  VB 
Ac"  — c  —  c  — c  —  c 


CAJ  *  L  -UR.  +  L  •  VB. 
Ad  — c  —  d  — c  —  d 


Enter  the  following  linear  correction  equation  into  the  system  of  such 
equations: 


Ap  =  CA  - -  +  CA  AaM  +  CA  Aa  +CAlI  AU 

An  n  A  a  „  xN  Aa  „  yH  AU  o 

-  o  xN  o  yH  *  o  o 

n  J 


+  CAo  Afl  +  C  A  ,  Ai  +  C  4  ...Ac"  +  C  .  .  tA 


dSra ijfetoJ&mficmfy 
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2.3.6  Altazimuth  Observations 

If  the  azimuth,  A,  and  elevation  angle,  h,  are  observed,  compute 
the  unit  vectors  S,  E,  and  Z: 


S  *  sin  0  cos  9 
x 


(a) 


sin  0  sin  9 


) 


S,  Southward  Unit  Vector 


S  *  -  cos  0 
z 


E  *  -sin  6  ' 

X  ! 


(b)  Ey  *  cos  9 


E,  Eastward  Unit  Vector 


E  -  0 
z 


Z  *  cos  0  cos  9  \ 


(c)  Zy  *  cos  0  sin  0  ^ 


Z,  Zenithal  Unit  Vector 


Z  *  sin  0 
z  ^ 


After  this,  compute  ,  A^  and  ^  ,  where  the  o  subscript  denotes  quantities 


o  o  o 

which  are  calculated  from  observed  data: 


»  -  cos  A  cos  h  \ 
o 


(d)  L  .  «  sin  A  cos  h  \ 

y  o  / 


L  ,  «  sin  h 

zh 

o 
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A  .  =  sin  A 

xn 

o 

A  ,  =  cos  A 

yho 

o 


D  ,  =  cos  A  sin  h 

xh 

o 


D  ,  =  -  sin  A  sin  h 

yh 
J  o 


D  ,  =  cos  h 

zh 

o 


Next  rotate  the  components  of  1^  >  ^  to  the  equatorial  coordinate 

system:  o  o  o 


(g)  L  =  L,S  +  L,E  +  L,Z 
&  -o  xh  -  yh  —  zh  — 
o  o  o 


A  =  A,S  +  A,E  +  A,Z 
— o  xh  —  yh  —  zh  — 

o  J  o  o 


D  =  DlS  +  D.E  +  D1Z 
— o  xh  —  yh  —  zh  — 

o  o  o 


Compute  AL  =  L  -  L 
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Form  the  partial  differential  coefficients  as  in  (2.3.5)  with  ^replacing 
L  and  enter  the  following  linear  correction  equation  into  the  system  of 
such  equations: 


Pc  <&-W-CAu  4^+  C  A*  +  CAVA>0  +  CA»A". 

no 


+  CAD  +  CA1  A1  +  CA."Ac''  +  CAdAd 


Again  form  the  coefficients  as  in  (2.3.5),  this  time  with  replacing  V 

and  enter  the  following  linear  correction  equation  into  the  system  of  such 
equations : 

p  (D  .  AL)  -  C  Ano  +  C.  Aa„  +  C.  Aa  -►  C  AlI  AU 

c  — o  An  -  Aa  „  xN  Aa  ..  yN  AU  o 

JL2.  n  xN  o  yN  o  o 

no 


+  cAfl  Afl  +  G  Ai  Ai  +  CAc„Ac»  +  CAdAd 


2.3.7  Equatorial  Angular  Coordinates 

If  a.  the  topocentric  right  ascension,  and  5  the  topocentric 
declination  are  observed,  compute  the  vectors  to  JL,  A,  and  D: 


L  =  cos  <5  cos  a 

X 


(a)  L„  *  cos  6  sin  a  ^  L 


L  =  sin  6 
z 
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D  * 
x 

o 


(C)  D  - 

Jo 


D 

z 

o 


-  sin  6  cos  a 

-  sin  <5  sin  a 

cos  6 


D 

— o 


Compute  aL  =  -  L^,,  then  form  the  partial  differential  coefficients  and 

compute  the  linear  correction  euqation  as  in  (2.3.6)  substituting  A  for 

and  D  for  D  . 

— o  — o 


2.3.8  Range-Rate  Observations 


(a) 


If  p,  the  slant  range  rate,  is  observed,  compute 
the  station  velocity  R  and  then  p  i 


where  6  =  0.058,834,47 


£c  -  i  +  4 
K  =  K  *  ic 
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2 

(b)  Compute  ex^  *  a  (e  cos  E-e  ) 

and  ey^  =  a  '•A.-e2  e  sin  E 

(c)  Compute  the  partial  differential  coefficients 
R±  and  U±: 

(1)  Ru->/f>a3/2  ssjjr'3 

(2)  R,  -  §  +  <0  -  VJ  Ry 

(3)  R^  =  ( *T ua5^2  r-3)  [sin  (E+(D)  -  a^  a  sin  E  -  a .  ] 

(4)  R  a  (*/Va5^2  r-3)  [-cos  (E+(n)  -  a  e  sin  E  +  a^  ] 

2R 

(5)  ic„_a 

3R 

(6)  Rd  -  -f  (U  -  u/ 

n 

(7)  Uu  =  -Jil  a3/2  ey^  r"3 

(»  i>n  -  f  +  (u  -  u0>  (io 

2 

(9)  a5/2  r'3)  /l-e2  [cos  (E+o>)  -  (1  +  ^  )  ] 

(IQ)'  Uyjj  =  WjU  a5^2  r  3)  ./  1-e2  [sin  (E+co)  -  (1  +  ^  )  ] 

2U 

(11)  U  =  —  (U  -  U  ) 

c  n  o 

3U  2 

(12)  ud  =  -y5-  (u  -  uor 

n 
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(d)  -Form  .the  gartialvcoefficicnfs  i 


C  .  =  L  •  Ur  P  (R.i  -  v  U  )  -  p  Rl  +  b  -  UR 

An  — c  —  L  c  n  n  rc  n  J  r-c  —  n 


L  •  V  fp  (U  +  -  U  )  -  p  U  1  +  p  •  VD 

— c  —  L  'c  v  n  r  n  'c  n  J  £-c  —  n 


C  Aa  =  -c  *  -  [  *>„<»«  -  4  V 

xN 


-  Pc  RxN  ]  +  £c  '  2  R: 


xN 


+  ic  •  X  [  <4*  +  f  <W  -  Pc  UXN  ]  +  K  •  2 


V  u 


xN 


CA3yN  "  ic  *  2  [  Pc  <RyN  -  ♦  V  -  4  RyN  3  +  ^  RyN 


c  v  yN  yN  rc  yN 

ic  *  ^  [  pc  ( V  +  f  V  '  4  V  1  +  4  *  x  V 


C  AlI  =  L  •  U  [p  (R  -  v  U  )  -  p  R  ]  +  P  •  DR, 
AU  —  c  —  L  *0  u  u  c  u J  —  u 

o 


+  L 
— c 


V  [  o  (u  +  -  U  )  -  p  U  ]  +  a  -  VU 
—  i-^c  u  ru  cuJ  c  —  u 


CAft  =  "  it*  £  pr  r  v  cos  i  +  L  •  V  cos  i  [  pc  r  -  Pc  r  1 


+  p  •  V  r  cos  i  +  L  •  W  sin  i  f  P  (r  v  sin  u  -  r  cos  u) 
c  —  — c  —  c 


p  r  cos  u  1  -  £  •  W  r  sin  i  cos  u 

c  J  c  — 
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^  ■  1^. '  W  [  Pc  (r  v  cos  u  +  i  sin  u)  -  Pc  r  sin  u  ] 


+  £  '  W  r  sin  u 

C  — 


C  Ac"  =  ic  •  £pc  <V*  V  ■  ^c  RJ  +  4 ' S  Rc 

+  Lc  •  V  [pc  (Uc  +  f  Uc)  -  pc  Uc]  +  £e  •  V  Uc 


Ad 


=  1^  -u[pc(Rd-vud)  -  PcRd]  +  4  *2 


^c  ‘  V-  [Pc  <”d  +  7  V  -  Pc  Ud  1  +  4  *  I  Ud 


(e)  Compute  A p  *  Pq  -  p  ,  the  slant  range- rate  residuals 

(£)  Enter  the  following  linear  correction  equation  into  the  system 
of  such  equations: 


p  Ap=  C 


An 
_ c 

An  n 


+  C 


Aa  „  +  C  A_  Aa  „  +  C  A„  AU 


AaxN““xNo 


%  X  AUo-° 


+  CA  a  Aflo  +  C At  Al  +  c  Ac"  Ac”  +  CAd  Ad 


2.3.9  Rejection  of  Observations 

In  the  normal  mode  of  operation,  this  program  rejects  observations 
in  the  same  manner  as  the  operational  differential  correction  program  (SGPDC).* 
The  absolute  maxima  for  residuals  of  range  and  angle  observations  are  set  to 
1000  km  and  for  range-rate  observations  to  0.5  km/sec  during  the  initial 
iterations  in  the  differential  correction  process.  The  relative  maxima  are 
computed  as  the  product  of  a  constant  (1.5)  times  the  computed  ms  errors  for 
the  combined  range,  angle  and  range-rate  observations.  (The  rms  is  computed 
as  the  square  root  of  the  sum  of  the  squares  of  the  accepted  residuals 
divided  by  the  square  root  of  the  number  of  accepted  residuals.)  When  the  rms 
of  the  positional  observations  become  less  than  50  km.,  the  absolute  rejection 
criterion  is  reduced  to  75  km. 


*  Aeronutronic  Publication  U-1691,  revised  1  October  1962. 
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In  the  override  mode,  where  a  specified  set  of  elements  is  corrected 
"n"  times  regardless  of  convergence  (see  fields  1  and  2  on  D.C.  control  card 
Figure  27),  the  observation  rejection  process  is  specified  on  the  D.C.  control 
card  (field  9)  in  one  of  two  ways: 

(1)  Absolute  maximum  rms  multiplier  method:  An  initial 
absolute  maximum  for  range  and  angle  observations  and 
an  rms  multiplier  for  each  correction  are  specified. 

(2)  Absolute  maximum  method:  An  absolute  maximum  for  range 
and  angle  observations  is  entered  for  each  correction 
specified. 

In  either  case,  the  rejection  of  range-rate  observations  is  treated  as  in  the 
operational  D.C.  program. 

A  further  modification  to  the  rejection  process  is  available  in 
both  the  normal  and  override  modes.  This  modification  is,  the  time  factor, 
t  -  t 

— -  ,  (where  t  -  tQ  is  in  units  of  days)  multiplying  all  absolute  maxima 

prior  to  the  rejection  test.  This  factor  is  not  used  if  it  is  less  than 
unity.  This  option  is  specified  in  field  seven  of  the  D.C.  control  card 
(Figure  27). 


2.3.10  Corrected  Elements 

Compute  the  corrected  elements  L  ,  a  „  ,  a  „  ,  h  ,  h  ,  h  ,c",d 
v  o’  xN’yN  x  y  z  ’ 

a  J  o  o  o  o 

N 

(a)  Let  Y,  C.  A.  =  v. ,  j  =  1,2,3, ..  .represent  all  of  the  linear 
i=l  13  1  3 

correction  equations  (i.e.,  the  C^'s  are  the  coefficients,  the  A^'s  are  the 
corrections  to  the  orbital  parameters  at  time  tQ,  the  v's  are  the  observation 
residuals,  and  N  is  the  number  of  parameters  being  corrected) 

Statistical  weighting  of  the  input  observations  is  accomplished  in 
this  program  by  assigning  a  constant  standard  deviation,  o  ,  for  each  quantity 
that  is  observed  by  a  particular  sensor.  These  standard  deviations  are  entered 
into  the  linear  correction  equations  is  shown  below.  They  are  introduced  to 
the  program  by  the  SIGMA  cards  (see  figure  14).  The  standard  deviation  for 
each  observation  may  be  varied  by  entering  four  multiplying  factors  on  each 
observation  card,  (see  figure  28) . 
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The  following  matrix  equation  is  solved  to  give  the  corrections, 
in  the  least  squares  sense,  to  the  orbital  parasMters  at  time  t  .  a,  is 
the  resumed  standard  deviation  of  the  observation  producing  v^.  ^ 


£  2  -2  y,  -2 

j  Clj°j  j  CljC2j0j  *  ' 


y  -2  y  2  -2 

j  CljC2j  °J  j  C2j  0j  *  * 


-2 


j  Clj°Nj  °5 


y  -2 

J  C2jSj 


^CC  n  If  (>  (J  /»2  Q  “2 

j  1 j  Nj  °j  j  C2jCNj  j  •  *  '  j  Slj  j 


y  _  -2 

j  cuvj  ai 


HjVj'2 


(b)  The  resulting  corrections  are  applied  as  follows  (a  prime  means 
that  the  element  is  a  corrected  element) ; 


An 

n'  -  n  (1  +  — —  ) 
o  o  n  ' 

o 


(c")'  =  c"  +  Ac" 


d'  =  d  +  Ad 


c.  .  .  LsDl  jd 

o  ,2  ’  360 

n 

o 


U'  =  U  +  AU 
0  0  o 
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xN  xN 
o  o 


+  A  a 


xN. 


"  ayN  +  ASyNc 
J  o 


Cl  '  =  SI  +  ACl 

o  o 


i  '  -  i  +  Ai 


L*  *  u'  +  n ' 

o  o  o 


If  «' 
z 


L' 

o 


U'  -  12' 
o  o 


if  V 

z 


=  a'2  +  a'2 

*N0  V 

k 


o  o 
\2/3 


a'  = 
o 


n1 


o  / 


i-f  J2 


Po  =  a0  (1  "  eo2> 


W  =  sin  Cl'  sin  i' 
xo 


W' 


cos  12  '  sin  i' 


W 


W'  =  cos  i' 
z 

o 


*  cos  i'  0 

*  cos  i'  <  0 


P 
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h'  -  fp*  W  \ 

xo  4  °  xo 

h'  =  (T7  W  ) 

yD  J  °  y0  f 


i'  «  p '  W  J 

z  ro  z  ' 

o  \  o 


W 

»•  = 

o  -W 


a 

_i  yN 

CD  *  s  tan  - — 

°  a^t 


where  the  quadrant  is  determined 
from  the  signs  of  the  numerator 
and  denominator 
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SECTION  3 

PROGRAM  OPERATION 


3.1  GENERAL  DESCRIPTION 

The  equations  appearing  in  Section  2  are  combined  in  a  Philco 
2000  computer  program.  The  program  is  divided  into  three  main  sections: 
a  differential  correction  subroutine,  an  ephemeris  calculation  sub¬ 
routine,  and  an  ephemeris  calculation  subroutine  provided  with  an  error 
analysis  section.  With  the  program  options  available,  these  sections 
can  be  used  in  any  combination  or  order. 

When  operating  this  program,  the  following  tapes  must  be  used 
as  shown  in  the  table,  regardless  of  the  section  of  the  program  being 
used. 


Table  III  Tape  Locations 


Tape 

Location  (Logical  tape  unit) 

Output  tape* 

5 

Standard  library  tape 

7 

Input  tape 

8 

SEAI  tape  (if  needed) 

10 

*N0TE:  Use  data  select  9  to  print  output  from  ephemeris 
section  of  program.  Use  data  select  0  to  print 
output  from  differential  correction  section  of 
program. 
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3.1.1  Program  Input 

As  input,  the  program  expects  the  N,  M  elements,  Lq,  a^N  ,  a  , 

a  7  °  7  ° 

and  h^;  the  drag  parameters  cq  and  d;  the  area  to  mass  ratio  “jjj-;  plus 

other  data,  depending  on  the  program  option  chosen.  This  is  discussed 
more  fully  in.  Section  3.3. 

3.1.2  Program  Output  (See  Section  3.3  for  samples  of  output) 


The  output  obtained  depends  upon  which  section  of  the  program 
is  being  used  and  also  upon  several  output  options. 


(a)  If  differential  correction  is  being  performed,  the  program  will 
print  out  the  correction  residuals,  the  old  elements  and  the  new  (corrected) 
elements,  along  with  many  auxiliary  quantities.  The  program  will  also 
punch  out  a  set  of  corrected  element  cards.  If  the  output  option  I0UT=1 
is  specified  with  differential  correction,  the  program  will  output  many 
quantities  calculated  in  the  ephemeris  subroutine.  These  include  the  com¬ 
ponents  of  r  and  r;  the  secular,  long-period,  short-period,  and  overall 
changes  in  the  elements;  the  values  of  the  89  general  perturbation  terms 
(see  section  3.6  for  units  of  terms)  and  also  many  intermediate  quantities. 


(b)  If  ephemeris  calculation  is  being  performed,  the  program  will 
output  rA  0^,  rA  ©3,|  A  £  |  >  Ar>  and  the  quantities  specified  by  the  out¬ 
put  option  I0PT  (see  summary  of  program  options,  Section  3.2)  after  each 

time  increment  from  t  to  t„._. 

o  END 

(c)  If  ephemeris  calculation  with  position-error  analysis  is  being  performed, 
the  program  compares  the  results  obtained  by  using  all  of  the  terms  with 

the  results  obtained  when  terms  selected  by  the  operator  are  omitted.  The 

program  then  outputs  the  following  quantities:  A  r,  rA  rA  A£  |, 

A  (2>  Au  and  A  i  (all  calculated  as  in  Section  2.2.8)  and  the  quantities 
specified  by  the  output  option  I0PT  after  each  time  increment  from  tQ  to 

^ND* 

3:2  PROGRAM  OPTIONS 

The  following  outline  is  a  summary  of  the  options  available  to 
the  usejr  of  this  computer  program.  The  application  of  these  options  is 
discussed  more  fully  in  Section  3.3. 
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3.2.1  I SENT 

(a)  ISENT  =  +1:  Input  sensor  data  with  cards  and  then  input 

SIGMA  cards  (Note:  This  must  be  +1  whenever  the 
D.  C.  Program  gets  its  sensor  data  from  cards,  even 
if  an  ephemeris  calculation  appears  first.) 

(b)  ISENT  =  -1:  Bypass  SIGMA  card  input  and: 

(1)  For  D.  C.  option,  input  sensor  data  from  tape 

(2)  For  ephemeris  option,  this  indicates  no  sensor  data  input 

(c)  ISENT  =  0:  (Blank  card)  Input  SIGMA  cards  and  input  sensor 

data  from  tape 


3.2.2  IDCEPH 

(a)  IDCEPH  =  0:  End  program  (Use  blank  card  for  0) 

(b)  IDCEPH  =  +1:  Go  to  differential  correction  subroutine 

(c)  IDCEPH  =  -1:  Go  to  ephemeris  calculation 

3.2.3  ICAL:  Used  Only  Under  Ephemeris  Option:  IDCEPH  =  -1 

(a)  ICAL  =1:  Go  to  ephemeris  calculation  without  error  analysis 

(b)  ICAL  =2:  Go  to  ephemeris  calculation  with  error  analysis 

3.2.4  I0PT:  Used  only  under  ephemeris  option:  IDCEPH  =  -1 


(a) 

I0PT  =  1: 

Output 

t, 

x,  y. 

z,  x,  y,  z 

(b) 

I0PT  =  2: 

Output 

t. 

Lat. , 

Long.  ,  H 

(c) 

I0PT  =  3: 

Output 

t. 

x,  y, 

z,  x,  y,  z  and  t,  Lat.,  Long.,  H 
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(d)  I0PT  =  4:  Output  t,  x, y,  z,  x,  y,  z  and  term  values  from  General 

Perturbations  section  of  XYZ  Subroutine 

(e)  I0PT  =  5:  Output  t,  Lat. ,  Long.,  H  and  term  values  from  General 

Perturbations  section  of  XYZ  Subroutine 

(f)  I0PT  =  6:  Output  t,  x,  y,  z,  x,  y,  z  and  t,  Lat.,  Long.,  H  and 

term  values  from  General  Perturbations  section  of  XYZ 
Subroutine 


3.2.5  IBACKl:  Used  Only  Under  Ephemeris  Option :  [IDCEPH  =  -1. 

(a)  IBACKl  =11:  Input  new  time  values  and  compute  another 

ephemeris  using  the  same  General  Perturbations 
terms  specified  •  by  NTEKMS(I)  (Note: 

If  IBACKl  =  1  is  used  after  the  Error  Analysis 
Section,  the  program  will  input  a  new  omitted 
term  case  and  compare  this  with  the  same  nominal 
term  case  used  with  the  preceding  error  analysis.) 

(b)  IBACKl  =  2:  Input  new  times  and  terms  and  compute  another 

ephemeris 

(c)  IBACKl  =  3:  Input  a  new  value  for  IDCEPH  to  determine 

whether  to  end  program  or  start  another  case 
with  new  elements 

(d)  IBACKl  =  4  :  End  program 


3.2.6  I0UT 

(a)  I0UT  =  0:  There  will  be  no  XYZ  output  from  D.C.  Subroutine, 

only  D.C.  output 

(b)  I0UT  =  1:  There  will  be  XYZ  output  from  D.C.  Subroutine  after 

each  batch  of  observations 


3.2.7  AG0M 

(a)  AG0M  ^  1:  Include  calculation  of  radiation  pressure  effects 

(b)  AG0M  <  1:  Omit  calculation  of  radiation  pressure  effects 
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3.2.8  DTERM 

(a)  DTERM  =  0:  Calculate  DTERM  in  the  BEGIN  subroutine 

(b)  DTERM  ^  0:  Use  this  value  for  d  and  bypass  the  calculation 

of  DTERM 


3.2.9  TERMS  (Two  cards:  Columns  1  through  80  on  the  first  card 

correspond  to  terms  Q ( 1)  through  Q(80)  in  the  program. 
Columns  1  through  9  on  the  second  card  correspond  to 
terms  Q(81)  through  Q(89)  in  the  program.) 

(a)  A  one  punched  in  any  of  the  above  mentioned  columns  will  cause 
the  corresponding  term  to  be  included  in  the  ephemeris  calcula¬ 
tion. 

(b)  A  zero  (or  no  punch)  punched  in  any  of  the  above  mentioned 
columns  will  cause  the  corresponding  term  to  be  set  equal  to 
zero  in  the  ephemeris  calculation. 

(c)  Identification  of  terms  Q(l)  through  Q(89)  is  found  in  Section  3.6. 


3.2.10  For  description  of  further  options  see  Differential  Correction 
Control  Card  format  (Figure  27) 


3.3  EXAMPLES  OF  PROGRAM  OPERATION 

The  data,  the  type  of  cards  contained  in  the  input  deck  and  also 
the  output  obtained  from  the  program  depend  upon  the  program  options  chosen. 
Therefore,  this  section  is  divided  into  three  subsections.  Each  subsection 
describes  one  sample  case. 

The  understanding  of  the  following  descriptions  will  be 
facilitated  by  referring  to  the  flow  chart  of  the  main  program  (Section  3.5). 

3.3.1  SAMPLE  CASE  #1:  Perform  differential  correction  of  elements 
with  sensor  data  input  from  cards. 
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(a)  Input  cards  used:  The  following  is  a  description  of  the  cards  used  and 
the  reasons  for  their  use,  listed  in  the  order  in  which  the  cards  appear 
in  the  input  deck.  (see  Figures  4  and  5) 

(1)  ISENT  CARD  (Figure  11)  with  a  1  punched  in  Column  6.  This 

card  tells  the  program  to  take  the  sensor  data  from  cards  (See  Section  3.2.1). 

(2)  3  SENS0R  CARDS:  These  cards  contain  sensor  data  in  standard 
SPADATS  format,  Figure  12.  The  three  cards  indicate  that  only  three  sensors 
are  used  for  this  case.  (Note:  no  more  than  200  sensor  cards  may  be  used) 

(3)  ENDSENS  CARD  (Figure  13) :  This  card  is  used  to  tell  the  pro¬ 
gram  that  all  SENS0R  CARDS  have  been  read  in. 

(4)  ENDSIGMA  CARD  (Figure  15) :  This  card  is  normally  used  after 
SIGMA  CARDS,  but  in  this  case  indicates  that  no  SIGMA  CARDS  are  used.  (Note: 
no  more  than  100  sigma  cards  may  be  used) 

(5)  IDCEPH  CARD  (Figure  16) :  With  a  1  punched  in  column  6.  This 
card  directs  the  program  to  the  differential  correction  subroutine  (see 
Section  3.2.2). 

(6)  7  ELEMENT  CARDS  (see  Figures  17  to  23) :  These  cards  contain 
the  initial  elements  needed  for  the  program.  They  are  in  standard  SPADATS 
format.  An  E  must  appear  in  column  80  of  each  card.  These  standard  elements 
are  converted  to  the  N,  M  elements  by  the  program. 

(7)  AG0M  CARD  (see  Figure  24  and  Section  3.2.7):  The  0.0  in 
columns  8,  9  and  10  is  the  DTERM  input.  Since  it  is  zero,  it  tells  the 
program  to  calculate  DTERM,  the  period  decay  acceleration. 

The  1  in  column  20  is  the  I0UT  option.  This  tells  the  program 
to  output  data  (see  Section  3.2.6)  from  the  XYZ  subroutine. 

A  7' 

The  0.5  in  columns  29  and  30  is  the  value  of  — *-  and  since  it  is 

in 

less  than  1,  it  causes  the  program  to  omit  the  calculation  of  radiation 
pressure  perturbations. 

(8)  2  TERMS  CARDS  (see  Figures  25  and  26):  These  cards  tell  the 
program  which  of  the  89  General  Perturbations  terms  to  include  in  the  calcu¬ 
lations.  In  this  case,  all  89  terms  are  included. 

(9)  D.  C.  C0NTR0L  CARD  (see  Figure  27):  The  1  in  column  15 
tells  the  program  to  output  all  residuals.  The  0  in  column  79  tells  the 
program  to  omit  use  of  0's. 
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FIGURE  4.  INPUT  CARD  DECK  FOR  SAMPLE  CASE  #1 
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(10)  44  OBSERVATION  CARDS  (see  Figure  28) :  These  cards  contain 
observation  data  In  standard  SFADATS  format.  The  44  cards  Indicate  that  the 
data  from  44  observations  are  to  be  included  in  the  differential  correction. 

(11)  END  CARD  (see  Figure  29) :  This  must  appear  after  each  set  of 
observations.  In  this  case,  there  is  only  one  set. 

(12)  ENDOBS  CARD  (see  Figure  30):  This  card  tells  the  program 

that  there  are  no  more  sets  of  OBSERVATION  CARDS  to  read  in. 

(13)  IDCEPH  CARD  with  0  punched  in  column  6.  This  card  causes 

the  program  to  end  after  the  ENDOBS  card  is  read  in,  which  is  after  the 

differential  correction  of  the  elements. 

(14)  END  DATA  (see  Figure  34):  This  card  is  used  to  tell  the  com¬ 
puter  that  all  the  data  for  this  job  has  been  read  in. 

(b.)  Output  from  Sample  case  #1:  Figure  6  shows  the  output  obtained  from 
sample  case  #1.  This  includes  all  of  the  quantities  described  in  3.1.2a 
above . 


0 

E 

i: 

[ 

L 


3.3.2  SAMPLE  CASE  #2:  Perform  ephemeris  calculation  without  error 
Analysis.  Output  all  quantities  possible. 

(a)  Input  cards  used:  The  following  is  a  description  of  the  cards 
used  and  the  reasons  for  their  use,  listed  in  the  order  in  which  the  cards 
appear  in  the  input  deck,  (see  Figure  7) 

(1)  ISENT  CARD  (see  Figure  11)  with  a  -1  punched  in  columns  5 
and  6.  This  card  tells  the  program  to  bypass  the  input  of  the  SENSOR  and 
SIGMA  cards  since  these  are  not  needed  for  ephemeris  calculation. 

(2)  IDCEPH  CARD  (see  Figure  16) :  with  a  -1  punched  in  columns 
5  and  6.  This  card  directs  the  program  to  the  ephemeris  calculation 
section. 


L 

L 
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(3)  7  ELEMENT  CARDS  (see  Figures  17  to  23) :  These  cards  contain 
the  initial  elements  needed  for  the  program.  They  are  standard  SPADATS 
format.  An  E,must  appear  in  column  80  for  each  card.  These  standard  elements 
are  converted  to  the  N,  M  elements  in  the  program. 

(4)  I0PT,  ICAL  CARD  (see  Figure  31):  The  6  in  column  6  is  the 
I0PT  indicator.  Since  it  is  6  it  will  cause  the  program  to  output  all 
quantities  possible  (see  Section  3.2.4).  The  1  in  column  12  is  the  ICAL 
indicator.  This  directs  the  program  to  the  ephemeris  calculation  without 
error  analysis  (see  Section  3.2.3). 

(5)  2TEKMS  CARDS  (see  Figures  25  and  26  and  Section  3.2.9):  These 
cards  tell  the  program  which  of  the  89  (general  Perturbations  terms  to  include 
in  the  ephemeris  calculation.  In  this  case,  all  89  terms  are  included. 

(6)  Input  for  XYZ  Subroutine  (see  Figure  32):  The  first  six  .t 
quantities  are  the  epoch  time  and  are  punched  as  shown  in  Figure  32. 

The  0.5  in  columns  38,  39,  and  40  in  AG0M.  This  is  the  quantity 
A7 

-jjj-  and  since  it  is  less  than  1,  it  causes  the  program  to  omit  the  calcula¬ 
tion  of  the  radiation  pressure  perturbations. 

The  0.0  in  columns  44,  45  and  46  is  DTERM.  Since  it  is  zero,  it 
tells  the  program  to  calculate  DTERM,  the  period  decay  acceleration. 

The  1000.0  in  columns  55  through  60  is  DELTAT.  This  is  the  time 
increment  to  be  used  in  the  ephemeris  calculation.  The  1000.0  in  columns 
67  through  72  is  TEND.  This  is  the  time  interval  to  be  covered  in  the 
ephemeris  calculation.  Since  TEND  is  equal  to  DELTAT,  the  program  will 
calculate  the  ephemeris  for  only  two  times,  tQ  and  t  +  1000  minutes. 

(7)  IBACK1  CARD  (see  Figure  33  and  Section  3.2,5)  with  a  4  punched 
in  column  6.  This  card  tells  the  program  to  end  after  the  ephemeris  cal¬ 
culation. 

(8)  END  DATA  CARD  (see  Figure  34):  An  END  DATA  CARD  (as  in 
Sample  Case  #1)  must  appear  after  the  last  card  in  the  data  deck. 

(b)  Output  from  Sample  case  #2 

The  program  prints  out  the  quantities  described  in  Section  3.1.2b 
as  shown  by  Figure  8.  Only  the  output  for  t  =  t  +  1000  minutes  is  shown. 
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This  is  the  output  which  appears  on  data  select  nine  for  each  observation  used  in  the  differential 
correction  section  of  the  program.  See  section  3.6  for  units  and  meaning  of  quantities  shown  above 

FIGURE  6  (cont) 


FIGURE  7.  INPUT  CARD  DATA  FOR  SAMPLE  CASE  NO 
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3.3.3  SAMPLE  CASE  #3:  Perform  Ephemeris  Calculation  With  Error  Analysis. 

Output  All  Quantities  Possible. 

(a)  Input  cards  used:  This  is  a  description  of  the  cards  used  and  the 
reasons  for  their  use,  listed  in  the  order  in  which  the  cards  appear  in  the 
input  deck  (see  Figure  9). 

(1)  ISENT  CARD  (see  Figure  11)  with  a  -1  punched  in  columns  5 
and  6.  This  card  tells  the  program  to  bypass  the  input  of  the  SENS0R  and 
SIGMA  cards,  since  these  are  not  needed  for  ephemeris  calculation. 

(2)  IDCEPH  CARD  (see  Figure  16)  with  a  -1  punched  in  columns  5  and 
and  6.  This  card  directs  the  program  to  the  ephemeris  calculation  section. 

(3)  7  ELEMENT  CARDS  (see  Figures  17  to  23) :  These  cards  contain 
the  initial  elements  needed  for  the  program.  They  are  standard  SPADATS  for¬ 
mat.  An  E  must  appear  in  column  80  of  each  card.  These  standard  elements 
are  converted  to  the  N,  M  elements  by  the  program. 

(4)  I0PT,  ICAL  CARD  (see  Figure  31) :  The  6  in  column  6  is  the 
I0PT  indicator.  Since  it  is  6,  it  will  cause  the  program  to  print  out  all 
output  possible  (see  Section  3.2.4).  The  2  in  column  12  is  the  ICAL  indicator. 
This  directs  the  program  to  the  ephemeris  calculation  with  error  analysis 

(see  Section  3.2.3). 

(5)  2  TERMS  CARDS  (see  Figures  25  and  26,  and  Section  3.2.9): 

These  cards  tell  the  program  which  of  the  89  General  Perturbations  terms  to 
use  in  the  nominal  case  for  the  error  analysis.  In  this  case  all  89  terms 
are  used. 


(6)  Input  for  XYZ  Subroutine  (see  Section  3.2):  The.  first  six 
quantities  are  the  epoch  time  and  are  punched  as  shown  in  Figure  32. 

The  0.5  in  columns  38,  39  and  40  is  AG0M.  This  is  the  quantity 

and  since  it  is  less  than  1,  it  causes  the  program  to  omit  the  calcu¬ 
lation  of  the  radiation  pressure  perturbations. 

The  0.0  in  columns  44,  45  and  46  is  DTERM.  Since  it  is  zero, 
it  tells  the  program  to  calculate  DTERM,  the  period  decay  acceleration. 
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The  1000.0  in  columns  67  through  72  is  TEND.  This  is  the  time 
interval  to  be  covered  in  the  ephemeris  calculation.  Since  TEND  is  equal 
to  DELTAT,  the  program  will  calculate  the  error  analysis  for  only  two  times, 
t  and  t  +1000  minutes  (Note:  no  more  than  300  time  increments  may  be 
used  in  ?he  error  analysis  section,  i.e.,  TEND  must  be  less  than  300.) 

DELTAT 

(7)  2  NTERMS  (I)  CARDS  (see  Figures  25  and  26,  and  Section  3.2.9): 
These  cards  tell  the  program  which  of  the  89  General  Perturbations  terms 

to  use  in  the  omitted  term  case  for  the  error  analysis.  In  this  case, 
only  terms  1,  6,  42,  49,  50  and  88  are  included. 

(8)  IBACKl  CARD  (see  Figure  33  and  Section  3.2.5)  with  a  4  punched 
in  column  6.  This  card  tells  the  program  to  end  after  the  error  analysis. 

(9)  END  DATA  CARD  (see  Figure  34)  An  END  DATA  CARD  (as  in  Sample 
Case  #1)  must  appear  after  the  last  card  in  the  data  deck. 

(b)  Output  from  Sample  Case  #3 

The  program  prints  out  the  quantities  described  in  Section  3.1.2c, 
as  shown  by  Figure  10.  Only  the  output  for  t  =  tQ  +  1000  minutes  is  shown. 


3.4  Input  Card  Format 

This  section  describes  the  input  cards  used  in  this  program. 
Examples  of  the  use  and  contents  of  these  cards  may  be  found  in  Section  3.3. 
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FIGURE  9.  INPUT  CARD  DATA  FOR  SAMPLE  CASE  NO. 


section  3.6  for  units  and  meaning  of  quantities  shown 

FIGURE  10  (Cont.) 


Field 

1 


Columns 
1  -  6 


2-12 


7  -  72 


uli  «•>••> 


Description* 


Input  specification  (right  adjusted  integer) 

'4-1 "  means  input  sensor  data  with  cards  and  then  input 
SIGMA  cards. 

"-1"  means  bypass  SIGMA  card  input  and  input  sensor  data 
from  the  SEAI  (sensor,  element,  acquisition,  infor¬ 
mation)  tape  for  D.C.  option  or  input  no  sensor  data 
for  ephemeris  option 

'0"  (or  Blank  card)  means  input  SIGMA  cards  and  input 
sensor  data  from  tape 

Not  used  as  this  time. 


*Note:  See  summary  of  program  options. 


FIGURE  11  ISENT  CARD 
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Field 

Columns 

Description 

1 

1-4 

Sensor  Number 

2 

5-11 

0°  (+N)  -  latitude  (decimal  1 

3 

12  -  19 

A°  (+W)  -  longitude  (  " 

4 

20  -  25 

H  (meters)  -  altitude  (  " 

5 

26 

Classification 
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27  -  30 

Sensor  Type 

7 

31  -  34 

Previous  Sensor  Number 

8 

35  -  36 

Number  within  Sensor  Complex 

9 

37  -  54 

Name 

10 

55  -  56 

Equipment  type 

11 

57 

Continent 

12 

58  -  59 

Country  or  State 

13 

60  -  72 

Comments 

14 

73  -  79 

Not  used 

15 

80 

Card  type  =  S  (Punch  0,  2) 

between  cols.  7-8) 

"  "  15  -  16) 

after  col.  25) 


FIGURE  12  SENSOR  CARD 
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Description 
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1  -  7 

ENDSENS  is  punched  in  columns  1 
This  card  must  appear  after  the 

through  7. 
last  SENSOR  card. 

2 

8-80 

Not  used  as  this  time. 

FIGURE  13  ENDSENS  CARD 
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Field 


Columns 


5  -  12 
13  -  20 
21  -  28 


Notes: 


Description 

Sensor  number  (right  adjusted  integer) 

Range  standard  deviation  (s.d.)  (meters) 

Azimuth  or  right  ascension  (s.d.) (deg.) 

Elevation  or  declination  (is. d.  ) (deg.) 
c^:  Range -rate (s.d.) (meters/sec.) 

Next  sensor  number 

Four  eight-column  fields  corresponding  to  fields 
2-5  but  used  with  sensor  indicated  in  field  6. 

1.  All  entries  are  floating  point  numbers. 

2.  See  Differential  Correction  Control  card  (Fig.  27) 
for  SIGMA  card  options. 


FIGURE  14  SIGMA  CARD 
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Field  Columns  Description 

1  1-8  ENDSIGMA  is  punched  in  columns  1  through  8.  This 

card  must  appear  after  the  last  SIGMA,  card  or  after 
the  ENDSENS  card  if  no  SIGMA  cards  are  used.  It  is 
not  used  when  ISENT  *  -1  (See  summary  of  program 
options.  Sec.  3.2) 

2  9-80  Not  used  at  this  time. 


FIGURE  15 
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Columns  Description 

1-6  Program  option  specification  (right  adjusted  integer) 

'H-l "  means  go  to  differential  correction  subroutine 
"-1 "  means  go  to  ephemeris  calculation  subroutine 
’0"  means  end  program 

When  using  the  D.C.  portion  of  the  program,  an 
IDCEPH  =  0  card  is  used  as  the  last  input  card 
in  order  to  end  the  program. 

2-13  7-80  Not  used  as  this  time. 


Field 
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FIGURE  16  IDCEPH  CARD 
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Field 


Field 

Columns 

Description 

1 

1  -  3 

Satellite  number  -  justified  right 

2 

4-6 

Element  set  number  -  justified  right 

3 

7 

Not  used 

4 

8 

Card  Number  (Card  #  =  1) 

5 

9-18 

Satellite  name  for  Element  File  Update 

6 

19  -  22 

Not  used 

7 

23  -  36 

Nq  -  Epoch  revolution 

8 

37  -  50 

e  -  Eccentricity 

9 

51  -  64 

i  -  Inclination  (degrees)  (i  =  0  may  not  be  used  in 
this  program) 

10 

65  -  79 

Not  used 

11 

80 

Card  type 

E  -  Nodal  Elements 

FIGURE  17  ELEMENT  CARD  1 
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Field 


Field 


Columns 


4-6 


9-22 
23  -  36 
37  -  50 
51-64 
65  -  79 


Description 


Satellite  number  -  justified  right 
Element  set  number  -  justified  right 
Not  used 

Card  number  (Card  #  =  3) 

P  -  Anomalistic  Period  at  Epoch  -  days/rev. 
a 

-  Right  ascension  of  ascending  node  -  degrees 

-  Argument  of  perigee  -  degrees 
qQ  -  Perigee -distance- earth  radii 
Not  used 


Card  type 


E  =  Nodal  Elements 


FIGURE  19  ELEMENT  CARD  3 
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Field 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


Columns 
1  -  3 

4  -  6 

7 

8 

9-22 
23  -  36 

37  -  50 

51  -  64 
65  -  79 
80 


Description 

Satellite  number  -  justified  right 
Element  set  number  -  justified  right 
Not  used 

Card  number  (Card  #  =  4) 

2 

cq  -  Rate  of  change  of  period  -  days/ (rev) 

0  -  Time  derivative  of  right  ascension  of  ascending 

node  -  degrees/day  (Not  used  by  program) 

-  Time  derivative  of  argument  of  perigee  - 
degrees/day  (Not  used  by  program) 

Not  used 

Not  used 

Card  type 

E  =  Nodal  Elements 


FIGURE  20  ELEMENT  CARD  4 
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Field 

Columns 

Description 

1 

1  - 

3 

Satellite  Number  -  justified  right 

2 

4  - 

6 

Element  set  number  -  justified  right 

3 

7 

Not  used 

4 

8 

Card  number  (Card  #  *  5) 

5 

9  - 

22 

d  -  decay  acceleration  (floating  point  entry) 

6 

23  - 

79 

Not  used 

7 

80 

Card  type 

E  =  Nodal  Elements 

FIGURE  21  ELEMENT  CARD  5 
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Field 

1 

2 

3 

4 

5 


Columns 
1  -  3 

4  -  6 

7 

8 

9-22 


6  23-36 

7  37-50 

8  51-79 

9  80 


Description 

Satellite  number  -  justified  right 
Element  set  number  -  justified  right 
Not  used 

Card  number  (Card  #  =  6) 

a  -  semi-axis  major  -  Earth  radii  -  (Not  used 
by  program) 

P  -  Nodal  period  -  days/rev. 

N 

c  -  rate  of  change  of  nodal  period  -  days/  (rev) 
Not  used 
Card  type 

E  =  Nodal  elements 
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FIGURE  22  ELEMENT  CARD  6 
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Field  Columns 

1  1-3 

2  4-6 

3  7 

4  8 

5  9-22 

6  23-36 

7  37-50 

8  51-79 

9  80 


FIGURE  22  ELEMENT  CARD  6 


Description 


Satellite  number  -  justified  right 
Element  set  number  -  justified  right 
Not  used 

Card  number  (Card  #  =  6) 

a  -  semi-axis  major  -  Earth  radii  -  (Not  used 
by  program) 

P  -  Nodal  period  -  days/rev. 

N  * 

c„  -  rate  of  change  of  nodal  period  -  days/ (rev)‘ 
N 

Not  used 
Card  type 

E  =  Nodal  elements 
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Field 


Col  mans 


Description 


1 

1 

- 

3 

Satellite  number 

2 

4 

- 

6 

Element  set  number 

3 

7 

Not  used 

4 

8 

Card  number  (Card  #  =  7) 

5 

9 

- 

22 

Not  used 

6 

23 

- 

29 

Initial  Revolution,  decimal  may  be  punched  in  column  29 

7 

30 

- 

36 

Fn  (Final  Revolution);  decimal  may  be  punched  in 

column  36 

8 

37 

- 

50 

Expiration  date  of  Bulletin,  in  format:  YM1DDHHMNSS.SS, 

decimal  punched  in  column  48 

9 

51 

- 

58 

RMS,  in  format  XXXXX.XX;  decimal  punched  in  column  56 

FIGURE  23  ELEMENT  CARD  7 
(sheet  1  of  2) 
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Field 

10 

11 


12 

13 


Columns 

59  -  66 
67 


68  -  79 
80 


Description 

Number  of  observations  used  in  obtaining  RMS 
I  STOP 

Blank  or  0  =  correct  the  inclination 
element 

1  =  do  not  correct  the  inclination 

2  =  do  not  correct  the  drag  parameter 
4  =  correct  time  equation  only 

Not  used 
Card  type 

E  =  Nodal  Elements 
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FIGURE  23  ELEMENT  CARD  7 
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Field 


Columns 

1-10 


11  -  20 


21  -  30 


31  -  80 


Description 

DTERM,  "d"  drag  coefficient  (1/min^)  (floating  point 
constant  or  0.0).  If  DTERM  *  0.0,  then  "d"  is 
calculated  in  the  program. 

I0UT,  Differential  Correction  (D.C.)  output  specification 
(right  adjusted  integer). 

'4-1"  tells  D.C.  subroutine  to  output  values  calculated 
in  the  XYZ  subroutine. 

,'0"  tells  D.C.  subroutine  to  bypass  output  of  values 
calculated  in  the  XYZ  subroutine. 

AGtfM,  ,  used  in  radiation  pressure  calculation 
(floating  point  constant).  If  AG4fef  ^  1  the  XYZ  subroutine 
calculates  radiation  pressure  effects.  If  AGOtf  <1  the 
XYZ  subroutine  bypasses  the  calculation  of  radiation 
pressure  effects. 

Not  used  at  this  time. 


FIGURE  24  AG6M  CARD 
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Field 

1 

2 

3 


Columns 

1 

2 

3 


Description* 
Fixed  point  1  or  blank 
Fixed  point  1  or  blank 
Fixed  point  1  or  blank 


80 


80 


Fixed  point  1  or  blank 


*Hote:  Columns  1  through  80  correspond  to  terms  Q(l)  through  Q(80)  in  the 
program.  A  one  punched  in  any  colunn  will  cause  the  corresponding 
term  to  be  included  in  the  ephemeris  calculation.  A  blank  column 
will  cause  the  corresponding  term  to  be  set  equal  to  zero  in  the 
ephemeris  calculation. 


FIGURE  25  TERMS  CARD  1 
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Description 


1-8 


10  -  11 


12  -  13 


Correction  specification.  "1"  =  Correct,  ’T)"  *  Do  not 
correct  the  corresponding  element.  Elements  are 
a,  axN,  a^,  U,  ft  ,  i,  c*1  and  d  in  that  order.  If 

field  not  used,  program  follows  nominal  SPADATS  sequence.* 

Number  of  times  to  correct  elements  specified  by 
Field  1,  regardless  of  convergence. 

Correction  pattern  specification.  1  of  10  choices  by 
number.  1  means  correct  with  nominal  procedure.  Others 
not  used  at  this  time. 

Convergence  criterion  override  -  minimum  X  change  in 
RMS  for  convergence  (integer  <  99  and  used  only  when  a 

correction  is  specified) 


♦Aeronutronic  publication  U-1691,  revised  1  October  1962,  p.  3-65 


FIGURE  27  DIFFERENTIAL  CORRECTION  CONTROL  CARD 
(sheet  1  of  2) 


-100- 


FleM 

Columns 

Description 

5 

14 

Punched  card  output  type  -  ''0”  or  blank  “  standard 
SPADATS  Format.  "1"  -  N  M  Format. 

6 

15 

Residual  output  designator.  '0"  or  blank  *  output 
first  and  last  pass  residuals.  *1"  *  Output  all 
residuals. 

7 

16 

Rejection  criterion  time  factor  designator.  'O"  or 
ii  A"  m  Do  not  used  (time)  factor.  ”1"  *  use 

3 

time  factor  in  residual  rejection. 

8 

17  -  18 

Absolute  maximum  for  range-rate  residuals.  Decimal 
point  assumed  between  col.  17  &  18. 

9 

19  -  78 

Observation  rejection  overrides  for  range  and  angle 
observations  only.  Ten  fields  of  six  columns  each  are 
provided  for  specifying  an  RMS  multiplier  and  absolute 
maximum  or  an  absolute  maximum  (only)  for  each  iteration 
specified  by  field  two  in  the  differential  correction 
process.  See  Section  2.3.9  for  further  discussion  of 
observation  rejection. 

10 

79 

If'TC'or  blank,  then  do  not  use  °'s  specified  on  SIQfA 
cards.  If  *1,"  then  use  o's  on  SIGMA  cards.  If  "2", 
then  use  o's  on  SIGMA  cards  and  mult,  factors  on 
obs.  cards. 

11 

80 

Not  used  at  this  time. 

FIGURE  27  DIFFERENTIAL  CORRECT! ON  CONTROL  CARD 
(sheet  2  of  2) 
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Columns  Description 

1  -  3  Satellite  number.  Column  1  contains  a  minus  sign  if  this  is 

a  classified  observation:  +  or  -  are  not  allowed 

*  o  o 

4  -  5  Equipment  Type 

6  -  9  Station  Number 

10  Accuracy  or  Signal  Strength 

11  -  15  Date 

16-24  Time  (Z) 

25  -  30  Elevation/declination.  Column  25  can  be  overpunched  +  or  -  . 

31  -  37  Azimuth/right  ascension.  Column  31  can  be  overpunebed  +  or  - 

38-44  Slant  range  (KM) 


*  A  minus  overpunch  in  col.  31  indicates  fields  7  and  8  are  declination  and 
right  ascension, respectively. 


FIGURE  28  OBSERVATION  CARD 
(sheet  1  of  4) 
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Field 


Columns 


Description 


10 

45  -  53 

11 

54 

12 

55  - 

57 

13 

58  - 

59 

14 

60  - 

61 

15 

62  - 

63 

16 

64 

-  65 

17 

66 

-  69 

18 

70 

19 

71 

-  72 

20 

73 

-  74 

21 

75 

-  76 

22 

77 

-  78 

23 

79 

24 

80 

Range  rate  (KM/sec)  with  Implied  decimal  point  between  ^ 
columns  46  and  47  or:  maximum  frequency  shift  (cycles/sec  ) 
with  implied  decimal  point  between  columns  52  and  53. 


Code  for  field  10 


} 

1 


At  observation  time* 

Maximum 

Minimum 

Time  interval 


0  or  A  indicates  range  rate  in  field  10. 

1  indicates  max.  freq.  shift  in  field  10. 


or,  if  col.  58  contains  a  -  punch. 
Brightness  then: 

cols.  55-57  *  radar  cross-section 
in  meters ^ 

cols.  59-63  =  frequency  (dec.pt. 

between  60  and  61) 

(NOTE:  Not  used  by  SPS) 


Date  or  line  number 
Message  number 
Equinox 

Multiplying  factor  for  standard  deviation  of  range 
Multiplying  factor  for  standard  deviation  of  range  rate. 
Multiplying  factor  for  standard  deviation  of  azimuth  or 
right  ascension. 

Multiplying  factor  for  standard  deviation  of  elevation 
or  declination. 

Switch  indicator  used  by  manual  system. 

Card  type  (code  type  =  Any  numeric  between  0-9)  identifies 
an  Observation  card.  0  =  Unknown,  1-9  coded  according  to 
the  Association  Status  as  determined  in  Report  Association. 


FIGURE  28  OBSERVATION  CARD 

(sheet  2  of  4) 
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COLUMN  10  (ACCURACY) 


Either  accuracy  or  signal  strength  may  be  indicated  in  column  10, 
coded  according  to  the  following: 


If  type,  in  columns  4  and  5,  is  31  or  greater,  column  10  contains 
signal  strength.  If  type  is  30  or  less,  column  10  contains 
accuracy. 


Code  Figure 

Accuracy 

Signal  Strength 

0 

Normal  observations  made  under 
fair  conditions. 

Signal  strength 
good,  reliable 
measurement . 

1 

Observations  slightly  under  par  due 
to  outside  interference  (e.g.  some 
clouds,  reduced  visibility). 

Signal  fair. 

2 

Observations  only  poor  due  to 
outside  interference. 

Signal  weak, 
results  poor. 

3 

Only  estimates  possible  (mal¬ 
function  of  instrument.  Too  short 
time  of  object  seeing). 

Signal 

questionable. 

4 

Doubtful  observations,  unable  to 
verify  either  object  or  instrument 
behavior.  Observations  should  be 
considered  only  as  tentative. 

COLUMNS  55  -  63  (CROSS  SECTION-FREQUENCY/MAGNITUDE) 


The  block  containing  columns  55  through  63  is  a  dual  purpose  block 
where  cross  section  and  frequency,  or  magnitude  and  time  interval  are  indicated. 

In  order  to  specify  cross  section  and  frequency,  a  minus  is  used  in  column  58. 

No  sign  is  used  in  column  58  when  this  block  contains  magnitude  and  time  interval. 

Cross  section,  given  in  square  meters,  is  listed  in  columns  55 
through  57.  To  indicate  less  than  one  square  meter  cross  section, 
use  appropriate  numbers  and  a  minus  in  column  55,  thus  in  effect, 
putting  a  decimal  point  before  column  55.  For  larger  values  where 
three  digits  would  not  be  sufficient,  use  a  plus  in  column  55  to 
represent  ton  times  the  indicated  value  (adding  a  zero  to  the  value 
listed) . 


FIGURE  28  OBSERVATION  CARD 
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following: 


Frequency  in  megacycles,  is  listed  in  columns  58  through  63  with 
the  decimal  point  understood  to  be  located  between  columns  60  and 
61.  In  rare  cases  it  might  be  desirable  to  increase  the  range  of 
frequency  given  either  side  of  the  decimal  point.  To  do  this,  use 
a  minus  in  column  63  to  move  the  point  one  place  to  the  left,  or 
a  plus  in  column  63  to  move  the  point  one  place  to  the  right. 


COLUMN  70  (EQUINOX) 

Column  70  contains  year  of  Equinox  as  specified  by  the 


0  =  year  of  date 

1  =  1900 

2  =  1925 

3  =  1950 

4  =  1975 

5  =  2000 

6  =  1850 

7  =  1855 

8  =  1875 

9  =  to  list  actual  year,  if  not  provided 

above,  list  last  two  digits  of  year 
in  columns  71  and  72  and  use  a  minus 
in  column  70  for  18  and  a  plus  in 
column  70  for  19.  Example:  Equinox 
of  1961  would  contain  '4-61 "  in  columns 
70,  71,  and  72. 
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Field 

Columns 

Description 

1 

1-3 

END  punched  in  columns  1,  2,  3.  This  card  must 
appear  after  each  batch  of  observation  cards. 

2 

4  -  80 

Not  used  at  this  time. 
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Field 


B 

B 

B 

E 

Field 

Columns 

Description 

1 

1  - 

6 

END0BS  punched  in  columns  1 
This  card  must  appear  after 
following  the  last  block  of 
cards . 

through  6. 
the  END  card 
observation 

2 

7  - 

80 

Not  used  at  this  time. 

FIGURE  30  END0BS  CARD 
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Field 

1 


Columns 

1-6 


7  -  12 


13  -  72 


Description 

I0PT,  XYZ  output  option  selector  (right  adjusted  integer) 
"1"  means  output  t,  x,  y,  z,  x,  y, 

"2"  means  output  t,  lat.,  long.,  height  (h) 

"3"  means  output  t,  x,  y,  z,  x,  y,  z,  and  lat.,  long.,  h 
"4"  means  output  t,  x,  y,  z,  x,  y,  z  and  term  values 
from  General  Perturbations  Section 
"5"  means  output  t,  lat.,  long.,  h  and  term  values  from 
General  Perturbations  Section 
"6"  means  output  t,  x,  y,  z,  x,  y,  z  and  lat.,  long.,  h 
and  term  values  from  General  Perturbations  Section 
ICAL, error  analysis  selector  (right  adjusted  integer) 

"1"  means  proceed  to  ephemeris  calculation  without 
error  analysis 

"2”  means  proceed  to  ephemeris  calculation  with  error 
analysis 

Not  used  at  this  time 
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mm 
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mm 
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1 

Field 

1 

Columns 

1  -  6 

IYEAR: 

2 

7-12 

M0NTH: 

3 

13  -  18 

IDAY: 

4 

19  -  24 

IH0UR: 

5 

25  -  30 

MINUTE: 

6 

31  -  36 

SEC0ND: 

7 

37  -  42 

AG0M, 

8 

43  - 

54 

9 

55  - 

66 

10 

67  - 

78 

11 

79  - 

80 

S 


Epoch 


Description 

right  adjusted  integer 
right  adjusted  integer 
right  adjusted  integer 
right  adjusted  integer 
right  adjusted  integer 
floating  point  entry 
A  7 

~jr~  ,  used  in  radiation  pressure  calculation. 

If  AG<JM  ]>  1  the  XYZ  subroutine  calculates 
radiation  pressure  effects.  If  AG0M  <  1  the 
XYZ  subroutine  bypasses  the  calculation  of 
radiation  pressure  effects 

DTERM,  "d"  drag  coefficient  (floating  point  constant  or 

0.0)  If  DTERM  =  0.0,  the  "d"  is  calculated  in  the 
program. 

DELTAT,  A  t,  the  time  increment  to  use  in  XYZ  sub¬ 
routine  (floating  point  entry) 

TEND,  time  interval  to  cover  in  XYZ  subroutine  (floating 
point  entry). 

Not  used  at  this  time 
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Field  Columns  Description* 

1  1-6  IBACK1,  ephemeris  exit  option  (right  adjusted  integer) 

"1"  means  input  new  time  values  and  compute 
another  ephemeris  using  the  same  terms 
specified  by  NTERMS(I) 

"2"  means  input  new  times  and  terms  and 
compute  another  ephemeris 

"3"  means  input  a  new  value  for  IDCEPH  to 
determine  whether  to  end  program  or  start 
another  case  with  new  elements 

"4"  means  end  program 

2-12  7-72  Not  used  at  this  time. 


*Note:  See  also  Section  3.2 
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Field 


AERONUTRONIC  DIVISION 


3.5  Flow  Charts 


The  following  pages  contain  the  flow  charts  of  the  main  program 
and  of  the  more  significant  subroutines. 

Standard  SPADATS  Subroutines  are  used  at  many  places  in  the  program. 
The  flow  charts  for  these  are  not  shown  here  but  a  list  of  their  names  and 
reference  page  numbers  (SPADATS  manual*)  is  shown  below. 


Subroutine 

Page 

Subroutine 

Page 

ARCTAN 

4-3 

PANT 

4-29 

GL0P 

4-29 

SEPSUB 

4-8 

INITEL 

4-59 

SNSGET 

4-65 

KL0K 

4-9 

TLC 

4-15 

0BSGET 

4-23 

XSRCH 

4-69 

0BSL0D 

4-27 

TABLE  IV  SPADATS  SUBROUTINES  USED  IN  PROGRAM 


*Aeronutronic  Publication  U-1691,  Revised  1  Obtober  1962 
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FLOW  CHART  OF  MAIN  PROGRAM 


I  SENT :  0. 


/  INPUT 

>  -/  SENSOR  DATA 
'  *1  WITH  CARDS 


r  INPUT  > 
k  SIGMA  CARDS J 


INPUT 

IDCEPH 


^  INPUT 
ELEMENT 
.  CARDS 


UNPACK 

ELEMENTS 


^DCEPH^W  ^ 


y  INPUT  >. 

/IMP  BEGIN.  >. 

f  INPUT  \  / 

/  DTERM,  I0UT,  \ 

J  BEGIN  JMP  BEGIN;. 

(  NTERMS  (I)  \  J 

\  AG0M  j 

NTHCN 

\  (89  TERMS)  /"PV 

r  INPUT  \ 
DC  CONTROL  CARD) 
AND  OBSERVATION  J 
V  CARDS  / 


fSET  I0UT  -0) 


I  SENT <  0 


f  INPUT 
[l0PT,  ICAL 


1  W  — ICAL> - >(2 

-''iCAL  =  l\  /ICAL  -2 


f  INPUT  ' 
SENS0R  DATA  WIT 
l  TAPE 


CALL  DC 


''NEXT  CARD  \YES  /ISXNO 
AN  END0BS  y*-<JDIV  ■  Oy-#l 
V  CARD  /  \  / 


RUN  OUT  ALL 
OBSERVATION 
CARDS 
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I 


FLOW  CHART  OF  MAIN  PROGRAM  (Continued) 


£ 


(input  X 
N  TERMS  (1)1 


#  iyear.mJSth,  IDAY,  4 
f IH0UR,  MINUTE, SECOND  1 
[  DELTAT, TEND, DTERM, 
agOm 


/JMP  BEGIN. BEGINS 
l  JMP  BEGIN. NTHCN 


Calculate 

r  NOMINAL  OR  OMITTED^ 
f  TERM  CASE.  I.E.,CALL\ 
XYZ  SUBROUTINE 


I0PT 
1  or  4 


CALL  \ 
OUTPUT  V 

SUBROUTINE  JT 
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FLOW  CHART  OF  MAIN  PROGRAM  (  continued  ) 
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FLOW  CHART  OF  EPHEMERIS  INITIALIZING  SUBROUTINE 


BEGIN 

(SUBROUTINE) 


XYZSBI 


ZEIS'S, 

AG0H> 

AG0MT 


SIN  i 

o 

COS  i 


SIN  Q 
o 

cos  Q 


II  - 


WY 

WZ 

SIN  I 

C0SI 

XINCL 

SIN0 

C0S0 

XN0DE0 

FSQ 

F0 

IOMEGA 

RTESQ 

RTE0SQ 

A 

RTA 

XN0 

DTERM 

C 

00 


JUMP 

TO 

ICYZSBIR 


XINCL 

0.0 


XINCLI 

DLTL 


RETURN 


XYZSBIR 


(L  ) 
o  o 

2 

ao 

1/P 

ao 

Ja 

XLSUN0X 

-  A0SQ 

—  PA0 

-*•  RTAUP 

v  up 

kl 

XYZKI 

k2 

-*•  XYZK2 

o 

•  • 

’-*■  XYZK30 

SiiMBi 


RETURN 


RETURN 


S'  DO  \ 
WE  WANT  v 
RAD.  PRESS. 
V  CALC.  / 


G0  T0 
XYZSBR 


- 

0MGAS0 

SIN  “so 

— 

SIN0M 

SIN  (Os 

-4* 

SIN0MS 

SIN  2CDso 

-*■ 

SIN20M 

SIN  3  O)so 

— 

SIN30M 

coso)  so 

-4 

C0S0M 

cosco  so 

—4 

C0S0MS 

COS  2(0 

so 

— 

C0S20M 

COS  3  CO 

so 

- 

C0S30M 

-TAN  CO 

so 

- 

TAN0M 

TAN  3  0)so 

- 

TAN30M 

|cos  iQ| 

AC0SI 
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FLOW  CHART  OF  EPHEMERIS  SUBROUTINE  (CONTINUED) 
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FLOW  CHART  OF  EPHEMERIS  SUBROUTINE  (CONTINUED) 


FLOW  CHART  OF  EPHEMERIS  SUBROUTINE  (CONTINUED) 
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-►RDELT3 

Ie  1 

— ♦•ABSDRB 

-*SU 

r. 

-♦  R 

k 

rk 

rd0t 

(-)k 

— ►  RVD0T 

i 

-*  XINCLI 

SIN  ^ 

-*•  SINU 

COS 

-►C0SU 

CMfi.L 
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z 

M 

— t-XMX  ! 

X 

M 

— *-XMY 

y 

M 

— ►XMZ 

z 

W 

-♦WX 

x 

W 

-♦WY 

y 

w 

_*-wz 

z 

u 

—  UX  i 

x 

u 

— ►UY 

y 

u 

-*-UZ 

z 

V 

-+VX 

x 

V 

-+VY 

y 

V 

-►VZ 

z 

X 

-♦X 

y 

~*Y 

Z 

— »z 

X 

— »-XD0T 

y 

— ►YD0T 

Z 

— *-ZD0T 

/IS  \ 
I0UT-  -1 


NO  ^ 
RETURN 
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FLOW  CHART  OF  EPHEMERIS  SUBROUTINE  (CONTINUED) 


£ 


OUTPUT  THE  FOLLOWING 
QUANTITIES  FOR  OFFLINE 
PRINTING: 

SATELLITE  NUMBER, 
SATELLITE  NAME,  ELEMENT 
SET  NUMBER,  T  ,  TIME 
of  epoch,  the  Individual 
VALUES  OF  THE  89  GENERAL 
PERTURBATION  TERMS, 

*«, '  H'  a-»sL’  H' 

w1'’"1'  n‘L’ 

AS)  ,  <D  ,LT,n„  A7T, 
AM,iL,i,  Ai,eL> 
x,  y,  z,  x,  y,  z,  , 

Au,  rfc.  At,  A  r » 


A(rv) ,  rk,  rAS^ 
r  A  e3>  I  £  | 


I 


^  RETURN^ 
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FLOW  CHART  OF  RADIATION  PRESSURE  SECTION 
OF  EPHEMERIS  SUBROUTINE 
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FLOW  CHART  OF  RADIATION  PRESSURE  SECTION 
OF  EPHEMERIS  SUBROUTINE  (CONTINUED) 


©H 


DILTL  (1)  -  DILTL  (2) 
DLTL 
e  -*•  E 


W* 

W> 

W. 


•wx 

■WY 

■WZ 


i  -*■  XINCLI 
Q  XN0DE 


THE  F0LL0WING  ARE 
CALCULATED 

DICPI  (IK0UNT) 
DISPI  (IK0UNT) 
DILWX  ( IK0UNT) 
DILWY  ( IK0UNT) 
DILTL  (IK0UNT) 


RETURN 


YES 


DICPI  (1)  -  DICPI  (2) 
-*DECPI 

DISPI  (1)  -  DISPI  (2) 
-►DESPI 

DILWX  (1)  -  DILWX  (2) 
-♦OELWX 

DILWY  (I)  -  DILWY  (2) 
-►DELWY 


h© 
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FLOW  CHART  OF  OUTPUT  SUBROUTINE 
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FLOW  CHART  OF  ARCTANGENT  SUBROUTINE 


-126- 


FLOW  CHART  OF  DIFFERENTIAL  CORRECTION  (D.C.)  SUBROUTINE 


ENTER  HERE  FROM 
MAIN  PROGRAM 


SET  ERROR  EXITS.  PICK 
UP  INITIAL  ADDRESS  FOR 
OBSERVATION  BUFFER  FROM 
THE  LOCATION  LPOBLOC. 
ENTER  AKLOK  ROUTINE  TO 
GET  CURRENT  TIME.  SET 
UP  PAGE  HEADING. 


ENTER  OBSGET  ROUTINE  TO  GET 
NEXT  OBSERVATION,  WITH  ASSOCIATED 
SENSOR,  IN  CORE  FORMAT. 


SET  UP  OBSERVATION  FLAG  TO 
INDICATE  WHICH  QUANTITIES  WERE 
OBSERVED.  ENTER  OBVEC  ROUTINE 
TO  PROCESS  THE  OBSERVATION  AND 
SAVE  IT  IN  OBSERVATION  BUFFER. 


IF  RESOP  0,  SET  THE 
OUTPUT  FLAG  TO  PRINT 
RESIDUALS,  CORRECTIONS, 
AND  CORRECTED  ELEMENTS 
FOR  EVERY  ITERATION. 
SAVE  INITIAL  ELEMENTS. 


ENTER  EDATE  ROUTINE  TO 
GET  EARLIEST  DATE  THAT 
OBSERVATIONS  ARE  DESIRED 
FOR  THIS  SATELLITE.  ENTER 
OBSLOD  ROUTINE  TO  OBTAIN 
ALL  OBS.  BEYOND  THIS 
DATE  FOR  THIS  SATELLITE 
WITH  THE  ASSOCIATED 
SENSORS  OR  TO  INITIALIZE 
THE  OBSGET  ROUTINE  IF 
OBS  ALREADY  INPUT  BY 
MAIN  PROGRAM 


SET  THE  MAXIMUM  FOR 
RANGE  AND  ANGLE 
RESIDUALS  (MAX)  TO 
ABSMX  AND  FOR  RANGE 
RATE  RESIDUALS  (MAX2) 
TO  0.5  KM/SEC.  SET 
TO  PICK  UP  FIRST 
OBSERVATION  FROM 
OBSERVATION  BUFFER. 


RESTORE  NEXT  OBSERVATION 
FROM  OBSERVATION  BUFFER 


ENTER  COEFFI  ROUTINE  TO  ACCUMULATE 
THE  SUM  OF  THE  SQUARES  OF  THE  GOOD 
RESIDUALS . 


COMPUTE  a„ 


ENTER 


CCCA  TO  INITIALIZE  THE 
ROUTINE  THAT  DETERMINES 
WHAT  ELEMENT  SET  TO 
CORRECT  AND  WHEN  DIVER¬ 
GENCE  OR  CONVERGENCE 
HAS  OCCURRED  (CCCC) . 

SET  ABSMX  TO  1000  KMS . 


COMPUTE  RMS1  =  ROOT  MEAN  SQUARE 
OF  THE  GOOD  RANGE  AND  ANGLE 
RESIDUALS,  RMS2  =  ROOT  MEAN  SQUARE 
OF  THE  GOOD  RANGE  RATE  RESIDUALS. 
SET  MAX  TO  m  '  RMS1  AND  MAX  2  TO 
m  •  RMS2  (m  =  CONSTANT) . 

INITIALIZE  LEAST  SQUARES  MATRIX. 
SET  TO  PICK  UP  FIRST  OBSERVATION 
FROM  OBSERVATION  BUFFER. 
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FLOW  CHART  OF  D.C.  SUBROUTINE  (CONTINUED) 


SET  SWITCH  1  TO  1 


RESTORE  THE  NEXT  OBSERVATION 
FROM  THE  OBSERVATION  BUFFER 


ENTER  COEFF2  ROUTINE  TO 
ACCUMULATE  THE  SUM  OF  THE 
SQUARES  OF  THE  GOOD  RESIDUALS 
AND  TO  ENTER  THE  COEFFICIENTS 
IN  THE  LINEAR  CORRECTION  EQUA¬ 
TIONS  INTO  THE  LEAST  SQUARES 
MATRIX. 


PRINT  OLD  AND  NEW  RMS’S,  LAST 
CORRECTIONS  TO  THE  ELEMENTS, 
NUMBER  OF  RESIDUALS  USED  (AND 
REJECTED)  COMPUTE  AND  OUTPUT 
THE  OLD  ELEMENTS  AT  THE  NODAL 
CROSSING  TIME  CLOSEST  TO  EPOCH. 
COMPUTE  AND  OUTPUT  THE  NEW 
ELEMENTS  AT  THE  NODAL  CROSSING 
TIME  CLOSEST  TO  EPOCH. 


IF  OUTPUT  FLAG  IS  SET,  OUTPUT 
RESIDUALS  AND  ASSOCIATED 
INFORMATION  FOR  THIS 
OBSERVATION. 


COMPUTE  RMS1  AND  RMS2.  SOLVE 
THE  LEAST  SQUARES  MATRIX  TO 
GET  (AS  REQUIRED) 


AU  ,AQ  ,  A  i  ,  Ac  AND 
o  o  o  o 

Ad.  SET  TO  ZERO  THi.  CORREC¬ 
TIONS  TO  THE  ELEMENTS  THAT 
WHERE  NOT  BEING  CORRECTED. 
OUTPUT  THESE  CORRECTIONS  IF 
THE  OUTPUT  FLAG  IS  SET. 


ENTER  CCCC  ROUTINE  TO  DETER¬ 
MINE  WHAT  TO  DO  NEXT 


DIVERGED 


/  SWITCH  \ 
[1  (ORIGINAL-] 
l  LY  SET  TO  j 

v  i)  y 


TINUE 

RECTING 


COMPUTE  AND  OUTPUT  NEW 
ELEMENTS  FOR  TIME  OF 
NODAL  CROSSING  INTO  THE  PRE¬ 
DICTION  REVOLUTION.  PUNCH  THE 
7  NEW  NODAL  ELEMENT  CARDS  FOR 
THIS  NODAL  CROSSING  TIME  WRITE 
EITHER  "AFTER  RECYCLE  ALL 
ELEMENTS  CORRECTED"  OR  "ALL 
ELEMENTS  CORRECTED  IN  FIRST 
ATTEMPT" 


EXIT  TO  MAIN 


SET  OUTPUT 
FLAG 


ENTER  CREKT  ROUTINE  TO  COMPUTE 
THE  CORRECTED  ELEMENTS  PRINT 
THESE  NEW  ELEMENTS. IF  THE  OUT¬ 
PUT  FLAG  IS  SET 


SET  SWITCH  1  TO  2 


H 
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FLOW  CHART  OF  D.C.  SUBROUTINE  (CONTINUED) 


WRITE:  "ELEMENT  CORRECTION  NOT  CONVERGING, 
SATELLITE  RUN  -TERMINATED,"  RESTORE  LAST 
GOOD  SET  OF  ELEMENTS.  WRITE  THE  OLD  AND 
NEW  RMS'S,  LAST  CORRECTIONS  TO  THE  ELE¬ 
MENTS,  NUMBER  OF  RESIDUALS  USED  AND  NUMBER 
OF  RESIDUALS  REJECTED,  ALL  FOR  LAST  GOOD 
ELEMENTS.  COMPUTE  AND  OUTPUT  THE  LAST 
GOOD  SET  OF  ELEMENTS  AT  THE  NEW  EPOCH. 


I 


PUNCH  THE  7  NODAL  ELEMENT  CARDS  FOR  THE 
LAST  GOOD  ELEMENTS  AT  THE  TIME  OF  NODAL 
CROSSING  INTO  THE  PREDICTION  REVOLUTION. 
OUTPUT  ORIGINAL  ELEMENTS  AT  ORIGINAL 
EPOCH  AND  AT  NEW  EPOCH.  WRITE  "LIMITED 
CONVERGENCE  ELEMENTS  NOT  REPLACED." 
WRITE  "ELEMENTS  SUCCESSFULLY  CORRECTED: 
(LIST  OF  CORRECTED  ELEMENTS) 

4  0  -►  IDIV 


+ _ 

EXIT  TO  MAIN 
PROGRAM 
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FLOW  CHART  OF  D.C.  SUBROUTINE  (CONTINUED) 


E 

E 


FLOW  CHART  OF  D.C.  SUBROUTINE  (CONTINUED) 
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FLOW  CHART  OF  D.C.  SUBROUTINE  (CONTINUED) 


^COEF^  2  ^ 


ENTER  XYZSB  ROUTINE  TO  COMPUTE 
r,  r  PLUS  MANY  INTERMEDIATE 
QUANTITIES.  COMPUTE  L-Lq  PLUS 
THE  COEFFICIENTS  R„,  R„,  RxN, 

RyN  >  Rc  ’Uu  »Hi  ,UxN  ’  ^yN  >  u  • 
COMPUTE  ()c,  R 


L 

— c 


COMPUTE  THE  COEFFICIENTS  C  An, 

II 

CAaxN  2AayN’  C AU’ 

cAfi»  cAi  *  CAC"*  <&d  SPENDING 

ON  WHICH  ELEMENTS  ARE  BEING 
CORRECTED.  COMPUTE  A p.  IF 

|  Ap|<MAX,  ADD  (A p)2  TO  RMS1 
SUMMATOR  AND  ENTER  THE  ABOVE  CO¬ 
EFFICIENTS  INTO  THE  LEAST 
SQUARES  MATRIX. 


COMPUTE  AL.  COMPUTE  THE 
COEFFICIENTS  CAn  CAfl  , 

Tp  axN 

cAayN»  ^U*  CA n>  cAi»  cAc"> 

C^j  DEPENDING  ON  WHICH  ELEMENTS 
ARE  BEING  CORRECTED,  FOR  A  OR 
a.  COMPUTE  THE  A  OR  *  RESIDUALj 
IF  IT  IS  <  MAX,  ADD  ITS  SQUARE 
TO  THE  RSM1  SUMMATION  AND  ENTER 
THE  ABOVE  COEFFICIENTS  INTO  THE 
LEAST  SQUARES  MATRIX 


COMPUTE  THE  COEFFICIENTS  C^n, 


CA, 


axN’ 


Q^U>  PA£2> 


CAi,  CAc..,  q^d  DEPENDING  ON 

WHICH  ELEMENTS  ARE  BEING 
CORRECTED  FOR  h  OR  5  .  COMPUTE 
THE  h  OR  5  RESIDUAL.  IF  IT  IS 
<  MAX,  ADD  ITS  SQUARE  TO  THE 
RMS1  SUMMATION  AND  ENTER  THE 
ABOVE  COEFFICIENTS  INTO  THE 
LEAST  SQUARES  MATRIX 


£ 


FLOW  CHART  OF  D.C.  SUBROUTINE  (CONTINUED) 


COMPUTE  p-,  ex^  ,  ey^  AND  THE 
COEFFICIENTS  Ru>,^n>  ^xNj  ^yN»  ^c> 
Rd>Uu,  Un>  Ujjjj,  Uyg  AHD  Ud,  Uc. 

COMPUTE  THE  COEFFICIENTS 

C^£  *  C^axjj  ’  CAa  N  ’  CAu  ’  CAfi  * 
n  J 

CAi  *  CAc"  ’  CAd 


DEPENDING  ON  WHICH  ELEMENTS  ARE 
BEING  CORRECTED,  COMPUTE  THE  RESIDUAL 
A p  .  IF|Ap|  <  MAX2,  ADD  (A p)2  TO 
THE  RMS2  SUMMATION  AND  ENTER  THE 
ABOVE  COEFFICIENTS  INTO  THE  LEAST 
SQUARES  MATRIX 
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FLOW  CHART  OF  D.C 


-135 


(CONTINUED) 


FLOW  CHART  OF  D.C.  SUBROUTINE  (CONTINUED) 
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FLOW  CHART  OF  D.C.  SUBROUTINE  (CONTINUED) 
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PLOW  CHART  OF  D.C.  SUBROUTINE  (CONTINUED) 


NO  ♦ 

1 

^CORSPEC  - 

I  t  0  —  DIVFL 

'  T  ; 

Jno 

RMS1  -►  OLDSUM 

ES  Q  RMS  2  -  0  'y 
NO 

IrJC2  ~  0LDSUM2 


(TEST  2>— 0LDSUM 

(TEST  3>—  0LDSUM2 

ABSMX2  -*  ABSMX 

D/0  —DIVFL 

SVELT  +  7  -*•  0RM5 
SVELT  +  8  -*  0KMS2 


NEXT  ELEMENT  SET  TO 
CORRECT  FROM  KTABLE 

BUFFER  — ►  KNTRL. 

(SETRPT)  — ►RPT. 

NUMBER  OF  ELTS  TO  BE 
CORRECTED  ■*  C0UWTL 


C0UNTL  -►A 
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FLOW  CHART  OF  D.C.  SUBROUTINE  (CONTINUED) 
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FLOW  CHART  OF  D.C 


SUBROUTINE  (CONTINUED) 


COMPUTE  g0  AND  ke£BQ.  ENTER 
XYZSB  ROUTINE  TO  COMPUTE 
ORBITAL  ELEMENTS.  COMPUTE 
h  .  SET  THESE  ELEMENTS  AS 
THE  NEW  EPOCH  ELEMENTS  INTO 
THE  ELEMENT  RECORD  IN 
MEMORY.  ENTER  BEGIN  ROUTINE 
TO  COMPUTE  MORE  INITIAL 
ELEMENTS.  ENTER  CEL0UT 
ROUTINE  TO  COMPUTE  ORBITAL 
ELEMENTS  AT  CLOSEST  NODAL 
CROSSING. 


^  CEL0UT  ^ 


COMPUTE  o\>,  Eg,  1^,,  Mq  AND 


N,  V  COMPUTE  go 
ke£so  ENTER  XYZSB  ROUTINE 
TO  COMPUTE  ORBITAL  ELEMENTS 
AT  NODAL  CROSSING.  COMPUTE 


AMD 


h,  Pa,  a), 
CROSSING. 


PN  AT  NODAL 
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AERONUTRONIC  DIVISION 


3.6  Program  Symbol  Definitions 

This  section  shows  many  of  the  symbolic  locations  used  in 
program  and  the  quantities  which  they  contain.  The  units  of  the  89 
Perturbations  terms,  Q(l)  to  Q(89) ,  are  also  shown. 


SYMBOL 

CONTENTS 

SYMBOL 

A 

a 

(earth  radii  (E.R. ) 

C0SI 

ABSDRB 

1  Ar  | 

(km) 

C0SII 

AE 

a 

e 

(e.r.); 

COSIL 

AG0M 

A 

m  ' 

2 

(  -22-  ) 
gm 

C0SILI 

A0 

a 

0 

(E.R.) 

COSILS 

AXN 

axN 

SL 

C0SI4 

AYN 

Xl 

C0S2IL 

AXNL 

X 

COSN0D 

AYNL 

X 

C0S0M 

AXN0 

axN 

o 

C0S20M 

AYN0 

X 

COS30M 

AXNS 

X 

C0SU 

AYNS 

X 

C0SUN 

C 

c" 

(1/min) 

C0S2U 

C0 

c 

0 

2 

(days /rev  ) 

C0S3U 

C0SE0 

cos  (e+cjd) 

C0S4U 

the 

General 

CONTENTS 
cos  i 

o 

cos  i 

cos  i  _ 
oL 

cos  i 

2. 

cos  X  _ 
oL 

4. 
cos  X 

o 

cos  2i  _ 
oL 

cos 

COS  CD 

so 

COS  2  CD 

so 

COS  3  CD 

so 

cos  u 

cos  u 
n 

cos  2u 
cos  3u 
cos  4u 
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SYMBOL 

CONTENTS 

SYMBOL 

CONTENTS 

DELIE* 

A  i 

(rad.) 

ESINE 

e  sin  E 

DEm 

n0  AM 

ESQ 

2 

e 

DELN0D 

Afi 

(rad.) 

FL08 

A 

DELN0E* 

Aft 

II 

HX0 

h 

V 

(E.R.)^ 

A 

o 

DELRE* 

Ar 

(km) 

HY0 

h 

y 

If 

DELTAT 

A  t 

(min) 

HZ0 

o 

h 

* r 

II 

4> 

o 

DELTI 

A  i 

(rad. ) 

NDECAY 

nD 

DELTR 

A  r 

(E.R.) 

0MEGA 

CD 

O 

(rad.) 

DELTRD 

A  r 

(E.R. /kg  1  min) 

0MGAS 

CD 

S 

M 

DELTRV 

A  (r(r) 

(E.R. /k^1  min) 

OMGAS0 

CD 

SO 

II 

DELTU 

Au 

(rad.) 

0MGDT 

d  CD 

dt 

(rad.  /min.) 

DELUE* 

Au 

It 

P 

p 

(E.R.) 

DTERM 

d 

(1/min^) 

P 

PL 

II 

E 

e 

PI 

T\ 

EC0SE 

e  cos  E 

P0 

Po 

(E.R.) 

ELSQ 

2 

6  L 

Q0 

qo 

II 

E0 

e 

0 

-.(1) 

ft  ^ 

(rad. /min . ) 

E01 

(E+CD) 

(rad. ) 

Q(2) 

«  2 

II 

E0SQ 

2 

e 

o 

Q(3) 

II 

*These  symbols 

are  used  in 

the  Error  Analysis 

Section  of  the 

program 
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E 

SYMBOL  CONTENTS 

SYMBOL 

CONTENTS 

|j  Q(4) 

,fi^(Rad.  /min.) 

Q(23) 

a  (dimensionless) 
XNg 

Q(S) 

Q.  ^ 

II 

Q(24) 

a  " 

xN9 

L  Q(6) 

“l 

II 

Q(25) 

a  " 

xN,_ 

r  Q(7) 

co2 

II 

Q(26) 

10 

a  " 

xN-, 

Q(8) 

“3 

II 

Q(27) 

11 

a  " 

yNl 

E  Q<9) 

II 

Q(28) 

a  11 

yN2 

Q(10) 

f 

“5 

II 

Q(29) 

a  M 

yN3 

L  Q(H) 

n6 

(rad.) 

Q(30) 

X 

j  Q(12) 

°7 

It 

Q(31) 

Q(13) 

°8 

II 

Q(88) 

a  ** 

yNe 

[  Q(14) 

fi9 

II 

Q(32) 

3yN7 

Q(15) 

j 

nio 

II 

Q(33) 

l  Q(16) 

3xNL 

(dimensionless) 

Q(34) 

a  ** 

yN9 

j  ‘  Q(17) 

3xN2 

II 

Q(35) 

a  11 

yNio 

Q(18) 

axN3 

II 

Q(36) 

3yNll 

J  Q( 19) 

axN. 

4 

II 

Q(37) 

3yN12 

Q(20) 

j 

axN5 

I* 

Q(38) 

a  M 

yN13 

L  Q(21) 

axN6 

It 

Q(39) 

(rad /min.) 

Q(22) 

axN? 

II 

Q(40) 

m2 

i 

I 

I 

I 
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l 

SYMBOL  CONTENTS 

SYMBOL 

CONTENTS 

[j  Q(41) 

M3 

(rad. /min) 

Q(60) 

Uj_  (rad.) 

Q(42) 

r 

*1 

fl 

Q(61) 

U6  " 

L  Q(43) 

*2 

tl 

Q(62) 

r..  (earth  radii 
(E.R. )) 

r  Q(44) 

W3 

If 

Q(63) 

r2 

Li 

Q(45) 

% 

II 

Q(64) 

r3 

[  Q(46) 

*5 

ff 

Q(65) 

ir1  (E.R. /k  ^min) 
i.  c 

Q(47) 

‘  L1 

(rad. ) 

Q(66) 

r2 

L  Q(48) 

L2 

tl 

Q(67) 

r3 

r  Q(49) 

L3 

II 

Q(68) 

r4 

L  Q(50) 

L4 

II 

Q(69) 

^5  " 

j'  Q(51) 

S 

II 

Q(70) 

^6  " 

Q(52) 

L6 

II 

Q(71) 

*7  ” 

L  Q(53) 

H 

II 

Q(72) 

r8 

r  Q(54) 

i2 

II 

Q(73) 

*9  " 

L  Q(55) 

13 

fl 

Q(74) 

j  Q(56) 

U1 

II 

Q(75) 

rv2 

Q(57) 

T  ' 

U2 

II 

Q(76) 

rv3 

L  Q(58) 

U3 

II 

Q(77) 

-4  " 

r  Q(59) 

L 

I 

I 

I 

1 

u4 

II 
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SYMBOL 

CONTENTS 

SYMBOL 

CONTENTS 

Q(79) 

*6 

(E.R.  /k"1  min.) 
c 

RTE0SQ 

^  1-e2 

Q 

Q(80) 

-7 

II 

RTESQ 

Q(89) 

*3 

It 

RTESQ 

Q(81) 

(rad.) 

RTP 

“CT  (E.R.)* 

*o 

Q(82) 

n  12 

II 

RVD0T 

rv  (E.  R.  /kg  min) 

Q(83) 

13 

It 

RVD0T 

(rv)k 

Q(84) 

n  14 

II 

SINE0 

sin  (B+  to) 

Q(85) 

** 

II 

SINI 

sin  i 

o 

Q(86) 

15 

II 

SINILI 

sin  i 

Q(87) 

>6 

«l 

SINILN 

SIN  i  _ 
oL 

R 

r 

(E.R.) 

SINILS 

2  n 
sin  i  _ 
oL 

R 

rk 

II 

SIN2IL 

sin  2i 

oL 

RDELTl 

rA  9  x 

II 

SINN 00 

sin  ft 

RDELT3 

rA93 

II  • 

SIN0M 

sin  u> 
so 

RD0T 

i 

(E.R.  /k*1  min.) 
e 

SIN20M 

sin  2  (d 

so 

RD0T 

*k 

it 

SIM30M 

sin  3  cd 

so 

RN 

r 

(E.R.) 

SINU 

sin  u 

n 

RN 

r. 

II 

SINUN 

sin  u 

kn 

n 

RTA 

Va 

(E.R.)* 
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SYMBOL 

CONTENTS 

SYMBOL 

CONTENTS 

SIN2U 

sin  2u 

WY 

W 

y 

SIN3U 

sin  3u 

WZ 

W 

z 

SIN4U 

sin  4u 

X 

x  (E.R. ) 

su 

u 

(rad.) 

Y 

y  " 

su 

UK 

If 

Z 

z 

TEND 

t  -  t 

0 

(min. ) 

XD0T 

x  (E.R./k  \nin 

THGR 

e 

(rad.) 

YD0T 

y 

gr 

THRG0 
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